WL-TR-92-2098 


AD-A263  497 


PULSE  MITIGATION  AND  HEAT  TRANSFER 
ENHANCEMENT  TECHNIQUES 

VOL  4  -  TRANSIENT  BEHAVIOR  OF  HEAT  PIPE  WITH 
THERMAL  ENERGY  STORAGE  UNDER  PULSE  HEAT  LOADS 


L. C.  Chow 

M. J.  Chang 

University  of  Kentucky 

Department  of  Mechanical  Engineering 

Lexington,  KY  40506-0046 


AUG  1992 

FINAL  REPORT  FOR  07/01/87  -  07/31/92 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  IS  UNLIMITED. 


AERO  PROPULSION  AND  POWER  DIRECTORATE 
WRIGHT  LABORATORY 
AIR  FORCE  MATERIEL  COMMAND 
WRIGHT'  PATTERSON  AFB  OH  45433-6563 


NOTICE 


When  Government  drawings,  specifications,  or  other  data  are  used  for 
any  purpose  other  than  in  connection  with  a  definitely  Government-related 
procurement,  the  United  States  Government  Incurs  no  responsibility  or  any 
obligation  whatsoever.  The  fact  that  the  government  may  have  formulated  or 
in  any  way  supplied  the  said  drawings,  specifications,  or  other  data,  is  not 
to  be  regarded  by  implication,  or  otherwise  in  any  manner  construed,  as 
licensing  the  holder,  or  any  other  person  or  corporation;  or  as  conveying 
any  rights  or  permission  to  manufacture,  use,  or  sell  any  patented  Invention 
that  may  in  any  way  be  related  thereto. 


This  report  is  releasable  to  the  National  Technical  Information  Service 
(NTIS) .  At  NTIS,  it  will  be  available  to  the  general  public,  including 
foreign  nations. 


This  technical  report  has  been  reviewed  and  is  approved  for  publica¬ 
tion. 


MICHAEL  J.  MORGJ 
Project  Engineer 


E.  BEAM  ,  Chief 
rmal  Technology  Section 


Dopnly  Chitii 

Acruspito  fOivcf  Dii/ision 
Aejo  Prcp^isiu.'!  &  Fowtr  Uirecloipte 


If  your  address  has  changed,  if  you  wish  to  be  removed  from  our  mailing 
list,  or  if  the  addressee  is  no  longer  employed  by  your  organization  please 
notify  WL/POOS  ,  WPAFB,  OH  45433-  6563  to  help  us  maintain  a  current 
mailing  list. 


Copies  of  this  report  should  not  be  returned  unless  return  is  required  by 
security  considerations,  contractual  obligations,  or  notice  on  a  specific 
document . 


REPORT  DOCUMENTATION  PAGE 

form  Approved 

OMB  No  0704  0)88 

PuOht  ffO-  f '  1  fi.i  J»“''  'it  ihi’, '  I  ,lle»  tii.^n  of  nf  is  •  sliTn.Hi'd  tc;  t  hi  lur  rPSEJonsC.  in.  tud<ci-.j  'h#»  fimj’  Jnr  r» 

Qdthetino  jra  in'!  1  ^  the  Oata  neodeO.  md  vompietir'}  and  rpv*r/».nq  the  lOHfciiOm^f  information  <..'mments  » 

crMieftiOi^  .t  ro  .ni  ludihy  su<j9t*siii>fn  tof  redn.  u'-}  this  Dutdpo  t-j  Washinqfon  Headdu.irters  ServKes.  PireLto/otr  fi'* 

Davis  Miqh  /.av.  -j  1/0*1.  Afhnqton,  v  A  43Q2  and  to  the  off*-  e  o*  Manaqerneni  arid  ifiidge!.  Paperwooc  Redurtir.n  Ff«.| 

-'•'A'A  t  'Osi/  j.  Ill  'Os  j  n. ,  d.*M  souk  es. 

din.-}  th's  tout  7»*n  esl’n' iir  ■'  m,'  '’her  aspect  of  this 

inloKnafinn  I  );>»T.itii'ii  .ind  H*‘IK  <  f..  I/’S  l.'ftefsoo 
•I  t  ;0. '04  0  tSrt).  Washington  i -C  PDSOj 

j  2.  REPORT  DATE 

3.  REPORT  TYPE  AN')  OATES  COVERED 

1 

^■1 4  rif  .1  WKimifi-bBSkfiLlVArMB! 

4.  TITLE  AND  SUBTITLE  puLSE  MITIGATION  AND  HEAT  TRANSFER  ^  funding  numbers 
ENHANCEMENT  TECHNIQUES 

VOL  4  -  TRANSIENT  BEHAVIOR  OF  HEAT  PIPE  WITH  C  F33615-87 

_ THBmii-fiMBRgY-aiORAgfi  UMDER  gULSE  HEAT  LOADS  PE  63218 

6.  AUTHOR(Jj ,  c .  Chow  PR  D812 

M.J.  Chang  TA  00 


C  F33615-87-C-2777 
PE  63218 
PR  D812 
TA  00 
HU  08 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  AOORESS(ES) 

University  of  Kentucky 

Department  of  Mechanical  Engineering 

Lexington,  KY  40506-0046 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

UK-ME-92-04 


WRIGHT  LABORATORY 
AIR  FORCE  MATERIEL  COMMAND 
WRIGHT  PATTERSON  AFB  OH  45433-6563 
WL/POOS,  Attni  MORGAN  513-2552922 


10.  SPONSORING /MONITORING 
AGENCY  REPORT  NUMBER 

WL-TR-92-2098 


i2a.Dj|f^gQ^^^V|iy^iLifx,fj«j^EH|jLEASE;  DISTRIBUTION  IS 
UNLIMITED . 


IZb.  DISTRIBUTION  CODE 


13.  ABSTRACT  (MaKimum 200 words)  ^  novel  design  of  a  high-temperature  axially  grooved  heat 
pipe  (HP)  which  incorporates  with  thermal  energy  storage  (TES)  to  mitigate  pulse 
heat  loads  was  presented.  Phase-change  material  (PCM)encapsulated  in  cylindrical 
containers  was  used  for  thermal  energy  storage.  The  transient  responses  of  the 
HP/TES  system  under  two  types  of  pulse  heat  loads  were  studied  numerically.  The 
first  type  is  pulse  heat  loads  applied  at  the  heat  pipe  evaporator;  the  second  type 
is  reversed-pulse  heat  loads  applied  at  the  condenser.  The  transient  response  of 
three  different  HP/TES  configurations  were  compared:  (1)  a  heat  pipe  with  a  large 
empty  cylinder  installed  in  the  vapor  core,  (2)  a  heat  pipe  with  a  large  PCM 
cylinder,  and  (3)  a  heat  pipe  with  six  small  PCM  cylinders.  It  was  found  that  the  | 
PCM  is  very  effective  in  mitigating  the  adverse  effect  of  pulse  heat  loads.  The  1 
six  small  PCM  cylinders  are  more  efficient  than  the  large  PCM  cylinder  in  relaxing  I 
the  heat  pipe  temperature  Increase  under  pulse  heat  loads. 


THERMAL  ENERGY  STORAGE,  PHASE-CHANGE 
MATERIAL.  PULSE  HEAT  LOAD 


17.  StCURlTr  CLASSIFICATION  18  SECURITY  CLASSIFICATION  19.  SECURITY  ClASSIFI  ZATION  20.  LIMITATION  OF  ABSTRACT 
OF  REPORT  OF  THIS  PAGE  OF  ABSTRACT 

UNCLASSIFIED  UNCLASSIFIED  UNCLASSIFIED  UL 


NSN  7Sa0  0'  280  5S00 


St.indarcl  Mjrm  (Rov  H 

*  '  v..  •  hv  •IN''*  •f'f  ■  '  <  'M 

.'  IM  Ul/ 


TABLE  OF  CONTENTS 


CHAPTER 

1 

2 

3 


4 


PAGE 


INTRODUCTION  1 

1.1  Statement  of  The  Problem  1 

1.2  Description  of  The  Heat  Pipe  3 

1.3  Literature  Review  6 

AN  nPRQVED  ALTERNATING- DIRECTION- IMPLICIT  METHOD  13 

2.1  Introduction  13 

2.2  Mathematical  Formulation  15 

2.3  Results  and  Discussions  26 

2.4  Concluding  Remark  36 

NUMERICAL  MODEL  38 

3.1  Heat  Conduction  Through  Pipe  Vail  and  Vick  41 

3.2  Melting  and  Solidification  of  PCM  47 

3.3  One- Dimensional  Vapor  Flow  Model  55 

3.4  Coupling  of  Vapor  Flow  Vith  Evaporation 

and  Condensation  60 

3.5  Liquid  Pressure  Drop  63 

TRANSIENT  BEHAVIOR  OF  HP/TES  SYSTEM 

UNDER  PULSE  HEAT  LOADS  APPLIED  AT  EVAPORATOR  64 

4.1  Limitations  of  the  Heat  Pipe  64 

4.2  Operation  Under  Normal  Conditions  68 

4.3  Operation  Near  HP/TES  Limitation  84 

TRANSIENT  BEHAVIOR  OF  HP/TES  SYSTEM 

UNDER  REVERSED  PULSED  HEAT  LOADS  APPLIED  AT  CONDENSER  93 

5.1  Description  of  The  Problem  93 

5.2  Analytical  Model  95 


Hi 


97 


5.3  Results  for  Heat  Pipes  Vith  Adiabatic  Section 

5.4  Results  for  Heat  Pipes  Without  Adiabatic  Section  119 

6  CONCLUSIONS  I33 

APPENDIX  Listing  of  The  Computer  Program  I37 

REFERENCES 


.  —  f .» 


IV 


LIST  OF  FIGURES 


FIGURE  PAGE 

1.1  Three  different  HP/TES  configurations  2 

1.2  Details  of  a  conventional  heat  pipe  4 

2.1  The  f- factor  modified  ADI  Method  22 

2.2  Coordinate  system:  parallelpiped  27 

2.3  Average  temperature  error  for  a  cube  with 

constant  wall  heat  flux,  t=  2  30 

2.4  Average  temperature  error  for  a  cube  with 

constant  wall  heat  flux,  t=  10  32 

2.5  Variation  of  average  temperature  error  with  f  factor 

for  a  cube  with  constant  wall  heat  flux,  t=  10  33 

2.6  Average  temperature  error  for  a  cube  with 

constant  wall  temperature,  r=  0.2  35 

2.7  Average  temperature  error  for  a  cube  with 

constant  wall  temperature,  t=  1.0  37 

3.1a  Nodal  map  of  the  heat  pipe  -  side  view  39 

3.1b  Nodal  map  of  the  heat  pipe  -  end  view  40 

3.2  Combined  specific  heat  including  latent  heat  effect  across 

the  melting  temperature  for  a  phase- change  material  48 

3.3  Enthalpy  over  a  2AT  interval  across  the  melting  temperature 

for  a  phase- change  material  50 

3.4  A  typical  location  of  the  fusion  front  on  the  nodal  map  51 

3.5  Temperature  distribution  during  melting  of  a 

one- dimensional  half- space  (Ste=  1)  56 

3.6  Location  of  the  fusion  front  during  melting  of  a 

one- dimensional  half- space  (Ste=  1)  57 

3.7  Coupling  of  vapor  flow  with  evaporation  and  condensation  62 

4.1  Operation  limits  for  a  grooved  heat  pipe  without  anything 

installed  in  the  vapor  core  66 

4.2  Operation  limits  for  a  grooved  heat  pipe  with  an  empty 

cylinder  mounted  in  the  vapor  core  67 


v 


70 


4.3  Transient  response  of  heat  pipes  with  a  sudden  increase  in 

0 

heat  input  from  q=  4.3  to  10  V/cm 


4.4  Definition  of  Q.  and  Q.  72 

J.I1  u 

4.5  The  ratio  between  and  for  the  same  operating 

conditions  depicted  in  Fig  4.3  74 

4.6a  Axial  variation  of  vapor  pressure  at  t=  100  s  for  the 

same  operating  conditions  depicted  in  Fig  4.3  75 

4.6b  Axial  variation  of  vapor  temperature  at  t=  100  s  for  the 

same  operating  conditions  depicted  in  Fig  4.3  76 

4.7  Comparison  of  results  from  Fig  4.3  and  those  obtained 

using  a  larger  time  step  and  wider  grid  spacings  78 

4.8  Transient  response  of  heat  pipes  with  pulsed  heat  loads 

q=  10  V/cm^  79 

4.9  Portion  of  PCM  melted  versus  time  for  pulsed  heat  loads 

q=  10  V/cra^  80 

4.10  Transient  response  of  heat  pipes  with  periodic  pulsed 

heat  loads  q=  10  V/cm^  82 

4.11  Portion  of  PCM  melted  versus  time  for  periodic  pulsed 

heat  loads  q=  10  V/cm^  83 

4.12  Transient  response  of  heat  pipes  with  a  sudden  increase 

2 

in  heat  inputs  from  q=  4.3  to  20  V/cm  85 

4.13  Axial  variation  of  liquid  and  vapor  pressure  at  t=  12s 

for  the  same  operating  conditions  depicted  in  Fig  4.12  87 

4.14  Variation  of  for  the  same  operating  conditions 

depicted  in  Fig  4.12  88 

•  • 

4.15  The  ratio  between  and  for  the  same  operating 

conditions  depicted  in  Fig  4.12  90 

5.1  Schematic  diagram  of  the  cooling  system  for  a  power 

generator  94 


VI 


5.2  Transient  response  of  heat  pipe  with  Cj^QQp=  1000  J/K 

under  a  reversed  heat  load  98 

5.3  Variations  of  heat  input  and  heat  output  of  heat  pipes 

with  1000  J/K  under  a  reversed  heat  load  100 

5.4  Transient  response  of  heat  pipe  with  C2^Qp=  10000  J/K 

under  a  reversed  heat  load  103 

5.5  Variations  of  heat  input  and  heat  output  of  heat  pipes 

with  C]^QQp=  10000  J/K  under  a  reversed  heat  load  104 

5.6  Comparison  of  the  transient  response  of  heat  pipe  without 

PCM  predicted  by  lumped  model  and  finite- difference  method  106 

5.7  Vapor  mass  flow  rate  of  heat  pipes  under  the  same 

operating  conditions  depicted  in  Fig  5.2  107 

5.8a  Axial  variation  of  vapor  pressure  for  the  same 

operating  conditions  depicted  in  Fig  5.2  109 

5.8b  Axial  variation  of  vapor  temperature  for  the  same 

operating  conditions  depicted  in  Fig  5.2  110 

5.9  Transient  response  of  heat  pipes  under  a  reversed- pulse 

heat  load  112 

5.10  Portion  of  PCM  melted  versus  time  for  heat  pipes  under  a 

reversed- pulse  heat  load  113 

5.11  Transient  response  of  heat  pipes  under  a  periodic 

reversed- pulse  heat  load  with  time  period  of  2000  s  115 

5.12  Portion  of  PCM  melted  versus  time  for  heat  pipes  under  a 
periodic  reversed- pulse  heat  load  with  time  period  of  2000  s  116 

5.13  Transient  response  of  heat  pipes  under  a  periodic 

reversed- pulse  heat  load  with  time  period  of  600  s  117 

5.14  Portion  of  PCM  melted  versus  time  for  heat  pipes  under  a 
periodic  reversed- pulse  heat  load  with  time  period  of  600  s  118 

5.15  A  heat  pipe  without  an  adiabatic  section  operates  under 

a  partially  reversed  heat  load  120 

5.16  Transient  response  of  heat  pipes  under  a  reversed  heat 

load  which  covers  757.  of  condenser  121 

5.17  Variations  of  heat  input  and  heat  output  of  heat  pipes 

under  a  reversed  heat  load  which  covers  757.  of  condenser  123 

5.18  Transient  response  of  heat  pipes  under  a  reversed  heat 


vii 


load  which  covers  507i  of  condenser 


124 


5.19 

Variations  of  heat  input  and  heat  output  of  heat  pipes 
under  a  reversed  heat  load  which  covers  507  of  co-'denser 

125 

5.20 

Transient  response  of  heat  pipes  un’<»r  a  reversed  heat 
load  which  covers  257  of  condenser 

127 

5.21 

Variations  of  hr  input  and  heat  output  of  heat  pipes 
under  a  reversed  neat  load  which  covers  257  of  condenser 

128 

5.22 

Vapor  mass  flow  rate  of  heat  pipes  under  a  reversed  heat 
load  which  covers  257  of  condenser 

129 

5.23 

Transient  response  of  heat  pipes  under  a  reversed- pulse 
heat  load  which  covers  257  of  condenser 

131 

5.24 

Portion  of  PCM  melted  versus  time  for  heat  pipes  under  a 
reversed- pulse  heat  load  which  covers  257  of  condenser 

132 

NQHENaATURE 


c  specific  heat 

C  heat  capacitance 

hydraulic  diameter  of  vapor  flow 

E  truncation  error  function,  defined  in  Eq.  (27) 
f  vapor  friction  coefficient  or  f- factor 

h  average  convection  heat  transfer  coefficient 

H  enthalpy 

I,J,K  number  of  nodal  points  in  x,  y  and  z  directions 
k  thermal  conductivity 

1^  length  of  condensation  region 

L  half-length  of  parallelepiped 

L’  dimensionless  half-length  of  parallelepiped 

characteristic  length 

A  axial  local  vapor  mass  flow  rate 
Ma  Mach  number 

P  vapor  pressure 

q  surface  heat  flux 

q  dimensionless  surface  heat  flux 

Q  total  heat  rate 

r,^,z  space  coordinates 

Re  axial  Reynolds  number  of  vapor  flow 

Re^  radial  Reynolds  number  of  vapor  flow 

Ste  Stefan  number  of  phase  change  material 


ix 


t  time 

T  temperature 

U,V  temperature  at  intermediate  time  steps 

U  mean  axial  vapor  v...'  city 

V  radial  vapor  velocity 

x,y,z  space  coordinates 
X,Y,Z  dimensionless  space  coordinates 
0  thermal  diffusivity 

2 

6  central  difference  operator 

t  average  temperature  error 

7  ratio  of  specific  heats 

A  stability  parameter,  defined  i'  Eq.  (8) 

f>  density 

a  Stef an- Boltzmann  constant 

r  dimensionless  time  or  Fourier  number 

0  dimensionless  temperature 

(  amplification  factor  of  truncation  error  function 

Superscript 
n  time  index 

Subscripts 

a  analytical  solution 

c  condenser 

e  evaporator 

g  power  generator 

hp  heat  pipe 

i,j,k  mesh  point  indices  in  x,  y  and  z  directions 


X 


1 


or  in  r,  0  and  z  directions 
liquid  state 


loop 

rev 

s 

w 

x,y,z 

0 

1,^3 


liquid  sodium  loop 

reverse 

solid  state 

wall  surface  of  parallelepiped 
indicate  x,  y  and  z  directions 
initial 

indicate  x,  y  and  z  directions 


xi 


CHAPTER  1 


INTRODUCTION 

1.1  Statement  of  The  Problem 

Future  space  missions  will  require  thermal  transport  devices  with 
the  ability  to  handle  transient  pulse  heat  loads  [1] ,  especially 
evaporator  loads  with  high  peak- to- average  power  ratios  and 
reversed- pulse  heat  loads  at  the  condenser.  Utilization  of  systems 
which  are  designed  to  handle  pulse  heat  loads  is  impractical  because  of 
the  large  masses  required.  Incorporation  of  thermal  energy  storage 
(TES)  into  heat  pipe  rejection  systems  can  be  a  promising  method  to 
mitigate  such  pul^e  heat  loads. 

Tilton  et  al.  [2]  and  El-6enk  et  al.  [3]  have  examined  the 
transient  response  of  a  heat  pipe  under  external  thermal  loading  at  the 
condenser.  However,  they  did  not  offer  any  suggestion  for  dealing  with 
these  pulse  heat  loads.  Some  configurations  to  reduce  the  dangers  of 
pulse  heat  loads  have  been  proposed  by  Beam  [4]  and  Sheffield  [5] ,  but  a 
detailed  analysis  of  HP/TES  mitigation  techniques  has  not  been 
attempted.  These  concepts  must  be  tested  and  understood  so  that  they 
may  be  integrated  successfully  into  an  overall  thermal  control  system 
design. 

The  surf ace- area- to  volume  ratio  of  the  TES  elements  is  an 
important  parameter  in  the  melting  and  solidification  process  of  the 
PCM.  In  this  research,  the  transient  responses  of  three  different 
HP/TES  configurations,  as  shown  in  Fig  1.1,  were  tested  and  compared 
under  a  variety  of  heat  load  conditions. 
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1 .2  Description  of  The  Heat  Pipe 

The  heat  pipe  is  an  innovative  device  capable  of  transferring  large 
quantities  of  heat  through  relatively  small  cross-sectional  areas  with 
very  small  temperature  differences.  Moreover,  a  heat  pipe  requires  no 
external  power  input  to  sustain  its  operation.  The  concept  behind  the 
heat  pipe  was  first  suggested  by  Gauler  [6]  in  1942.  However,  it  was 
not  widely  publicized  until  1964  when  Grover  [7,8]  and  his  colleagues  at 
the  Los  Alamos  Scientific  Laboratory  independently  reinvented  the 
concept.  Grover  also  demonstrated  the  device  effectiveness  as  a 
high-performance  heat  transmission  mechanism,  coined  the  name  "heat 
pipe,"  and  developed  several  applications  for  such  systems.  Since  then, 
the  remarkable  properties  of  the  heat  pipe  have  become  appreciated,  and 
serious  developmental  work  is  still  taking  place. 

A  heat  pipe  consists  of  a  closed  tube  or  chamber  with  various 
shapes  whose  inner  surface  is  lined  with  a  porous  capillary  wick  as 
shown  in  Fig  1.2.  Wire  screen,  fiber  glass,  porous  metal,  and  woven 
cloth  have  all  been  used  as  the  capillary  wick.  Narrow  grooves,  cut 
lengthwise  in  the  interior  pipe  wall,  can  also  serve  as  a  capillary  wick 
structure.  The  wick  is  saturated  with  the  liquid  phase  of  a  working 
fluid  and  the  remaining  volume  of  the  tube  contains  the  working  fluid 
vapor  phase.  Heat  applied  at  the  evaporator  by  an  external  source 
vaporizes  the  working  fluid  in  that  section.  The  resulting  difference 
in  pressure  drives  vapor  from  the  evaporator  to  the  condenser,  where  it 
condenses  releasing  the  latent  heat  of  vaporization  to  a  heat  sink  in 
that  section  of  the  pipe.  Depletion  of  liquid  by  evaporation  causes  the 
liquid- vapor  interface  in  the  evaporator  to  enter  into  the  wick  surface 
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Heat  Input  Heat  Removal 


Evaporator  Adiabatic  Section  Condenser 


End  View 


Figure  1.2  Details  of  a  conveutional  heat  pipe 
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auld  a  capillary  pressure  is  developed  there.  This  capillary  pressure 
pumps  the  condensed  liquid  back  to  the  evaporator  for  reevaporation. 
Thus,  the  heat  pipe  can  continuously  transport  the  latent  heat  of 
vaporization  from  the  evaporator  section  to  the  condenser  section 
without  drying  out  the  wick.  This  process  will  continue  as  long  as  the 
flow  passage  for  the  working  fluid  is  not  blocked  and  a  sufficient 
capillary  pressure  is  maintained.  Heat  pipes  may  be  operated  over  a 
broad  range  of  temperatures  by  choosing  an  appropriate  working  fluid. 

The  amount  of  heat  that  can  be  transported  by  the  latent  heat  of 
vaporization  is  usually  several  orders  of  magnitude  larger  than  that 
which  can  be  transported  as  sensible  heat  in  a  conventional  convective 
system.  The  heat  pipe  can,  therefore,  transport  a  large  amount  of  heat 
with  only  a  small  unit  size.  Because  only  a  very  small  temperature  drop 
occurs  in  the  vapor  flow,  heat  pipes  have  thermal  characteristics  orders 
of  magnitude  better  than  any  known  solid. 

However,  imlike  solid  conductors,  the  heat  pipe  characteristics  are 
dependent  not  only  upon  size,  shape,  and  material  but  also  upon 
construction,  working  fluid,  and  heat  transfer  rate.  Moreover,  heat 
pipes  operate  under  heat  transfer  limitations  such  as  the  sonic  limit, 
capillary  limit,  entrainment  limit,  and  boiling  limit.  Vhen  any  of 
these  limitations  is  encoimtered,  the  capillary  wick  structure  may  dry 
out,  leading  to  failure  of  the  heat  pipe. 

In  this  research,  a  stainless  steel  grooved  heat  pipe  using  sodium 
as  the  working  fluid  was  studied.  Grooved  heat  pipes,  which  use  narrow 
grooves  as  the  capillary  wick  structure,  have  been  used  quite 
successfully  in  practice.  The  advantages  of  a  grooved  heat  pipe  are  its 
low  capillary  liquid  flow  resistance,  high  reliability,  and  ease  of 
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fabrication.  However,  the  rather  large  nominal  diameter  of  the  grooves 
makes  both  the  capillary  limit  and  entrainment  limit  quite  low.  Vinz 
and  Busse  [9]  and  Barantsevich  et  al.  [10]  showed  experimentally  that  if 
the  grooves  in  the  evaporator  section  are  wrapped  with  layers  of  fine 
wick  mesh,  the  capillary  limit  and  entrainment  limit  can  be  improved 
significantly  while  maintaining  low  capillary  liquid  flow  resistance. 


1.3  Literature  Review 

Analysis  of  transient  heat  pipe  behavior  has  been  the  subject  of 
many  recent  studies  [11-16].  Chang  and  Colwell’s  [11]  model  for  low 
temperature  heat  pipes  neglects  the  hydrodynamics  of  the  vapor  flow  and, 
therefore,  cannot  predict  the  vapor  pressure  and  temperature  variations 
along  the  heat  pipe.  Kuramae  [12]  and  Tilton  et  al.  [13]  did  not 
include  vapor  flow  in  their  models.  Bowman  et  al.  [14]  and  Jang  et  al. 
[15]  applied  too  complicated  vapor  flow  models  which  lead  to  lengthy 
computation  time  and  only  represent  the  earliest  stages  of  heat  pipe 
transient  behavior  effectively.  The  model  developed  by  Seo  et  al.  [16] 
is  impressive.  However,  this  two-dimensional  (axial  and  radial)  model 
cannot  predict  the  transient  behavior  of  a  heat  pipe  with  angularly 
nonsymmetric  boundary  conditions.  It  would  be  desirable  to  assimilate 
all  of  the  current  works  by  incorporating  the  strongest  features  of  each 
previous  approach  into  a  state-of-the-art  model. 

A  transient  three-dimensional  numerical  method  is  desired  to  model 
the  heat  conduction  through  the  heat  pipe  wall  and  wicks,  including  the 
liquid  in  the  grooves.  Thibault  [17]  gave  a  comparison  of  nine 
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numerical  schemes  for  the  solution  of  the  transient  three-dimensional 
heat  diffusion  equation.  The  major  drawback  of  the  pure  explicit  scheme 
is  that  the  stability  criterion  At/ (AX)  is  very  small  and  should  not 
exceed  0.5.  The  computations  may  become  prohibitively  expensive  because 
of  the  need  to  employ  very  small  time  steps.  On  the  other  hand,  the 
pure  implicit  scheme  cannot  be  used  efficiently  for  the  solution  of 
multidimensional  problems.  In  three  dimensions  it  necessitates  the 
solution  of  large  sparse  matrices,  which  requires  extremely  long 
computation  time.  After  considering  the  relative  accuracy,  computation 
time,  and  the  computer  core  storage  requirement,  Thibault  concluded  that 
alternating- direction- implicit  (ADI)  finite- difference  methods  offered 
the  best  compromise.  The  ADI  methods  only  require  solving  tridiagonal 
matrices,  so  they  have  a  tremendous  advantage  in  computation  time 
compared  to  the  pure  implicit  method.  Chang  and  Colwell  [11]  and  Jang 
et  al.  .[15]  used  two-dimensional  ADI  methods  to  model  heat  pipe 
transients  with  great  success.  The  conventional,  two-dimensional  ADI 
method  is  unconditionally  stable.  However,  when  extended  to  three 
dimensions,  the  conventional  ADI  method  is  conditionally  stable,  and 
very  small  time  steps  are  required  to  ensure  convergence  and  stability. 
In  heat  pipe  modeling,  a  small  Ar  is  needed  due  to  the  slender  geometry 
of  the  heat  pipe,  and  only  a  very  small  time  step  At  (about  0.001s)  can 
be  used  with  the  conventional  three-dimensional  ADI  method.  Chang  et 
al.  [18]  developed  a  new  three-dimensional  ADI  method  by  modifying  the 
conventional  ADI  method  with  an  f- factor  (0<f<l).  This  modification 
allows  the  time  step  to  be  increased  by  about  2  orders  of  magnitude 
without  significantly  compromising  the  accuracy  of  the  numerical 
solution.  They  also  showed  that  this  new  ADI  method  yields  much  higher 
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accuracy  than  the  well-known  Douglas’  [19]  and  Brian’s  [20]  ADI  methods. 
Details  of  this  new  ADI  method,  which  was  used  to  model  the  heat 
conduction  through  the  heat  pipe  wall  and  wicks,  will  be  introduced  in 
Chapter  2. 

In  the  numerical  solution  of  heat  conduction  problems  with  phase 
change  (Stefan’s  problem)  by  finite  differences,  enthalpy  methods 
[21-25]  or  heat  capacity  methods  [26-30]  can  be  used.  The  former 
methods  require  either  an  explicit  procedure  which  may  lead  to 
convergence  problems,  or  iteration  at  each  time  step  if  an  implicit 
procedure  is  used.  The  latter  methods  are  subject  to  the  problem  of 
jumping  the  latent  heat  peak,  necessitating  the  use  of  very  small  time 
steps  to  avoid  underprediction  of  the  phase  change  time.  Recently, 

Hsiao  [31,32]  proposed  a  new  finite- difference  method  for  Stefan’s 
problem.  In  his  scheme,  the  equivalent  heat  capacity  at  a  node  is  a 
function  of  the  temperature  at  that  node  and  all  the  surrounding  nodes. 
Hsiao  concluded  that  his  method  can  avoid  the  problem  of  jumping  the 
latent  heat  peak  and  allows  the  use  of  a  relatively  large  time  step. 
Hsiao’s  method  was  tested,  but  a  large  energy  balance  error  was  found. 
Pham  [33,34,35]  suggested  a  simple  method  which  includes  features  from 
both  the  enthalpy  and  heat  capacity  methods.  Comparing  this  method  with 
other  existing  methods  for  test  problems  with  exact  solutions,  Pham 
pointed  out  that  most  of  the  methods  agree  to  within  0.27.  with  the 
analytical  result,  except  for  Hsiao’s  method  which  yielded  results  with 
up  to  227.  error.  The  low  accuracy  of  Hsiao’s  method  could  be  due  to  its 
ambiguous  theoretical  basis.  Pham  also  concluded  that  his  method  is 
much  faster  than  other  methods.  However,  Pham’s  method  has  a 
singularity  problem  in  finding  the  equivalent  specific  heat.  In  this 
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research,  we  adopted  the  best  features  of  Phan’s  method  and  made  some 
modifications  to  improve  on  its  weak  points.  Compared  with  analytical 
solutions,  the  present  method  for  melting  and  solidification  was  found 
to  have  very  good  accuracy  without  the  singularity  problem  of  Pham’s 
method. 

The  vapor  flow  in  the  evaporator  of  a  heat  pipe  is  dynamically 
similar  to  pipe  flow  with  blowing  through  a  porous  wall,  while  the 
condenser  section  is  analogous  to  flow  with  suction.  The  axial  vapor 
mass  flow  increases  along  the  length  of  the  evaporator  region  to  a 
maximum  value  at  the  end  of  the  evaporator;  it  will  then  decrease  along 
the  condenser  region.  The  behavior  of  the  vapor  flow  in  a  heat  pipe  is 
also  similar  to  that  of  a  gas  flowing  through  a  converging- diverging 
nozzle.  In  the  heat  pipe  the  area  remains  constant  but  the  mass  flow 
varies.  In  a  nozzle,  the  mass  flow  is  constant  but  the  cross-sectional 
area  is  changed.  There  are  many  theoretical  investigations  of  heat  pipe 
vapor  dynamics  [36- 45] .  Most  of  this  research  work  focuses  on  the  vapor 
flow  itself  instead  of  studying  the  vapor  flow  phenomena  coupled  with 
heat  pipe  transients.  The  majority  of  these  vapor  flow  models  are  very 
complicated  and  impractical  because  their  calculations  require 
prohibitively  large  amount  of  computer  time.  Bowman  [46]  reached  a  very 
important  conclusion  that  since  the  response  time  of  the  vapor  dynamics 
is  very  short  compared  to  the  heat  transfer  response  time  of  the  heat 
pipe  wall  and  wicks,  the  vapor  flow  can  be  modeled  as  a  quasi- steady 
process.  He  also  studied  the  compressible  vapor  dynamics  in  a  heat  pipe 
experimentally  using  air  flowing  in  a  porous  tube  with  blowing  and 
suction  along  the  wall.  He  modeled  the  vapor  flow  using  a  numerical 
solution  to  the  axisymmetric,  unsteady  Navier- Stokes  equations  and  a 
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steady  one- dimensional  solution  technique.  Bowman  concluded  that  it  is 
adequate  to  treat  the  vapor  flow  as  steady  and  one- dimensional ,  and  he 
also  provided  suitable  vap^r  pressure  drop  correlations.  In  our 
research,  this  steady,  one- dimensional  vapor  flow  model  was  used  and 
coupled  with  the  evaporation  and  condensation  in  the  heat  pipe. 

The  transient  response  of  a  heat  pipe  under  normal  operating 
conditions  is  mainly  controlled  by  its  thermal  capacity,  conductance, 
and  vapor  temperature  drop,  and  is  only  slightly  influenced  by  the 
liquid  dynamics.  However,  liquid  dynamics  become  very  significant  when 
dryout  and  rewetting  occur  in  the  wick.  The  dryout  phenomena  can  cause 
a  dramatic  temperature  increase  at  the  evaporator  section  and  may  affect 
the  overall  heat  transport  device.  To  predict  dryout  of  the  wick,  a 
detailed  liquid  flow  model  in  the  wick  is  needed.  Beam  [47]  and  Ambrose 
[48]  et  al.  investigated  the  transient  behavior  of  the  liquid  flow  in  a 
heat  pipe  wick  using  a  one- dimensional  flat- front  model.  Their  models 
for  the  liquid  flow  utilized  a  simplified,  lumped- parameter  solution  to 
the  energy  equation  to  predict  the  temperature  response  and  determine 
the  mass  flux  of  vapor  out  of  the  evaporator  region.  Both  Beam  [47]  and 
Chang  and  Colwell  [11]  concluded  that  the  transient  response  of  the 
working  fluid  in  screen  wicks  is  so  fast  that  acceleration  terms  in  the 
momentum  equation  for  the  liquid  are  negligible.  However,  dryout 
depends  on  the  heated  zone  and  the  instantaneous  local  saturation. 

Thus,  if  dryout  is  to  be  accurately  predicted,  the  temporal  dependence 
of  the  saturation  distribution  must  be  taken  into  account.  Ambrose  [49] 
developed  a  technique  utilizing  X-ray  radiography  to  measure  liquid 
distribution  in  the  porous  wick  structures  of  a  heat  pipe  with  beryllium 
walls.  He  also  presented  a  new  transient  liquid  flow  model  to  predict 
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the  continuous  saturation  distribution  in  the  wick  structure.  The 
measured  saturation  distributions  compared  favorably  with  those 
predicted  by  the  liquid  model.  However,  solution  of  this  new  transient 
liquid  flow  model  requires  knowledge  of  the  saturation  dependence  of  the 
capillary  flow  properties,  which  can  only  be  determined  by  experiment. 
Kamotani  [50]  and  Hwangbo  et  al.  [51]  developed  two  very  similar  liquid 
flow  models  for  a  grooved  heat  pipe.  Their  models  can  predict  the 
variations  of  liquid  meniscus  contact  angle  and  liquid  pressure 
variation  along  the  heat  pipe.  However,  they  assumed  the  liquid 
meniscus  is  always  attached  to  the  top  of  the  groove  side  walls. 
According  to  the  experimental  verifications  given  by  Ogushi  et  al .  [52] , 
this  assumption  is  not  correct.  They  also  ignored  the  effect  of  vapor 
pressure  variation  on  the  liquid  meniscus  contact  angle.  This  is  not  a 
good  simplification  for  high  temperature  heat  pipes  where  the  vapor 
pressure  drop  is  always  the  dominant  one.  A  complete  capillary  liquid 
flow  model  for  a  screen  covered  groove  like  the  one  we  recommended  in 
this  research  is  urgently  needed  in  the  near  future  to  predict  the 
dryout  and  rewetting  behaviors  of  a  heat  pipe  with  this  type  of  wick.  A 
good  capillary  liquid  flow  model  should  include  the  effect  of  vapor 
pressure  changes  on  the  liquid  meniscus  contact  angle  and  also  be  able 
to  predict  the  saturation  distribution  analytically. 

Kamotani  [53]  analyzed  the  thermal  behavior  of  the  condenser 
section  of  a  heat  pipe  with  axial  rectangular  grooves.  Some  of  the 
vapor  condensation  occurs  on  the  land  areas  between  grooves.  The  liquid 
forms  a  thin  film  on  the  land  surface,  and  heat  is  removed  from  the 
vapor  through  the  liquid  film.  Vhen  the  heat  pipe  is  operating,  the 
liquid  in  each  land  area  is  drawn  continuously  into  the  grooves  by 
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capillary  force  and  then  transported  to  the  evaporator  section  along  the 
grooves.  Because  of  viscous  drag,  liquid  flow  in  the  grooves  suffers  a 
pressure  drop,  and  consequently  the  curvature  of  the  liquid  free  surface 
varies  along  the  path.  This  meniscus  variation  will  influence  the 
motion  of  the  liquid  film  on  the  land  and  thus  the  condensation  rate. 
Kamotani  claimed  that  the  thickness  of  the  liquid  film  on  the  land 
surface  is  only  on  the  order  of  several  microns.  In  this  research,  we 
neglected  the  thermal  resistance  of  the  condensed  thin  liquid  film 
because  it  is  miniscule  when  compared  with  the  resistance  of  the  heat 
pipe  wall. 
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CHAPTER  2 


AN  IMPROVED  ALTERNATING- DIRECTION- IMPLICIT  METHOD 


The  conventional  three-dimensional  alternating- direction  implicit 
(ADI)  method  is  modified  by  introducing  an  f- factor  (0<f<l).  This 
modification  allows  the  time- step  limit  be  increased  by  a  factor  of  1/f 
with  the  solutions  remaining  stable  and  retaining  high  accuracy.  This 
new  method  is  tested  for  two  different  boundary  conditions:  a  constant 
heat  flux  and  a  sudden  heating  of  the  surface  to  a  constant  temperature. 
In  addition,  it  is  compared  with  the  popular  Brian  and  Douglas  ADI 
methods.  Results  show  that  the  new  ADI  method  has  higher  accuracy  and 
requires  less  computer  storage  than  the  Brian  and  Douglas  ADI  methods. 


2.1  Introduction 

The  diffusion  of  heat  in  solids  has  numerous  applications  in 
various  branches  of  science  and  engineering.  Generally,  there  are  two 
different  approaches  to  deal  with  this  type  of  problems:  analytical  and 
numerical  approaches.  The  analytical  methods  are  usually  only 
applicable  to  linear  problems  with  simple  geometries.  On  the  contrary, 
the  numerical  methods  are  useful  for  handling  practical  problems 
involving  nonlinearities,  complex  geometries  and/or  complicated  boundary 
conditions. 

Thibault  [17]  gave  a  comparison  of  nine  numerical  schemes  for  the 
solution  of  the  transient  three-dimensional  heat  diffusion  equation. 
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Considering  the  relative  accuracy,  the  computation  time,  and  the 
computer  core  storage  requirement,  he  recommended  that 
alternating- direction  implicit  (ADI)  finite- difference  methods  are  among 
the  most  preferred  methods.  The  conventional  two-dimensional  ADI  method 
was  introduced  by  Peaceman  and  Rachford  [54]  in  1955.  The  advantage  of 
the  ADI  method  is  that  only  tridiagonal  matrices  need  to  be  solved. 
However,  when  extended  to  three  dimensions,  the  conventional  ADI  method 
is  conditionally  stable  and  very  small  time  steps  are  required  to  ensure 
convergence  and  stability.  Other  forms  of  the  ADI  method  include  the 
Douglas  method  [19]  and  the  Brian  method  [20] .  Douglas  and  Brian  ADI 
methods  are  unconditionally  stable,  they  possess  the  advantages  of  the 
implicit  scheme  with  no  limitation  on  the  size  of  the  time  step. 

However,  Thibault  [17]  pointed  out  that  these  two  unconditionally  stable 
ADI  methods  cannot  retain  good  accuracy  if  the  time  step  is  more  than  2 
times  larger  than  the  time- step  limit  required  for  the  conventional  ADI 
method. 

In  this  research,  the  conventional  three-dimensional  ADI  method  is 
modified  by  an  f- factor  (0<f<l).  A  very  important  characteristic  of 
this  modification  is  that  it  is  consistent  with  physical  considerations, 
and  is  not  just  based  on  mathematical  manipulations.  This  modification 
allows  the  time- step  limit  to  be  increased  by  approximately  a  factor  of 
1/f  without  compromising  significantly  on  the  accuracy  of  the  numerical 
solution.  This  new  ADI  method  is  presented  and  compared  with  the  Brian 
and  Douglas  ADI  methods  for  two  cases  where  analytical  solutions  are 
available.  Compared  with  the  Brian  and  Douglas  methods,  this  new  ADI 
method  has  higher  accuracy  when  large  time  steps  are  used.  Also,  the 
present  method  requires  less  computer  storage. 


14 


r 

2.2  Mathematical  Formulation 

First,  we  will  look  at  of  the  formulations  of  the  existing 
three-dimensional  ADI  methods:  conventional,  Brian  and  Douglas  methods. 
Then,  the  proposed  new  ADI  method  designed  to  overcome  the  shortcomings 
of  these  existing  ADI  methods  will  be  introduced. 

The  differential  equations  for  the  three-dimensional  heat  diffusion 
equation  can  be  written  as 


1  5T  _  5^T 

a  "dt  a.  2  ^ 


ah 


Introducing  dimensionless  parameters: 


X= 


,  Y=  ,  Z- 


z 


9= 


(1) 


Eq.  (1)  becomes 


99  _  ^9  ^9  ^ 


(2) 


Conventional  ADI  Method 

In  the  conventional  ADI  method,  the  heat  diffusion  equation  is 
solved  implicitly  in  turn  in  the  three  coordinate  directions  for 
one- third  of  the  tine  increment  each  [55].  The  basic  finite- difference 
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equations  for  each  of  the  three  one- third  time  steps  can  be  expressed 
as: 


.  -  +  6^^  .  ,  +  .  , 
ir/S  X  1,3.  y  i,j,k  z  i,j,l 


(3) 


=  sh.  .  ,  +  sh.  .  ,  +  sH.  .  , 

IrfZ  X  1,3 ,k  y  1,3 ,k  z  1,3 ,k 


(4) 


^+1  .y 

i>J  _  ^2y  ^  r2y  ^ 

St/S  X  i,3,k  y  i,3,k  z  i,3,k 


(5) 


For  convenience  of  analysis,  we  let  AX=AY=AZ.  After  rearranging 
Eq.  (3),  it  becomes: 


=  '^.j.k  ^  ^.j-i,k  *  «?.j.i.k  *  ^.j,k-i  *  »;,j,k.i  (6) 


Similar  equations  can  be  easily  derived  from  Eqs.  (4)  and  (5)  for 
the  y-  and  z-  directions.  Physically,  an  increase  in  the  central  nodal 
temperature  or  an  increase  in  any  one  of  the  neighboring  nodal 
temperatures  at  the  old  time  step  should,  with  other  conditions 
remaining  unchanged,  lead  to  an  increase  in  the  central  nodal 
temperature  at  the  next  1/3  time  step.  This  implies  that  all  the 
coefficients  on  the  right-hand  side  of  Eq.  (6)  must  have  the  same  sign 
(positive)  as  the  coefficient  of  U-  •  .  .  In  other  words,  negative 
coefficients  on  the  right-hand  side  of  Eq.  (6)  make  the  equations 
physically  unrealistic  and  may  lead  to  low  accuracy  [56] .  Same 
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statements  can  be  made  regarding  to  the  equations  for  the  y-  and  z- 
directions. 


On  the  right-hand  side  of  Eq.  (6),  only  the  coefficient  for  •  . 

could  be  negative  if  the  time  step  At  is  large.  In  order  to  have 
positive  coefficient  for  ,  it  is  required  that: 

<  0.75  (7) 

(AX)2 

Since  AX=AY=AZ,  the  equations  for  the  y-  and  z-  directions  require 
the  same  condition  in  Eq.  (7)  to  hold.  The  other  important  thing  we 
need  to  consider  is  the  stability  problem.  Ve  define  the  stability 
parameter  A  as  follows: 


A  =  (8) 

The  stability  criterion  for  the  conventional  three-dimensional  ADI 
method  is  [55] : 

A  <  1.5  (9) 

Equation  (7)  is  the  condition  for  the  solution  of  the  conventional 
three-dimensional  ADI  method  to  be  physically  realistic  and  have  good 
accuracy.  Equation  (9)  is  the  criterion  for  the  solution  to  be  stable. 
The  main  disadvantage  of  the  conventional  ADI  method  is  that  it  is 
conditionally  stable  and  very  small  time  step  is  required. 
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Brian’s  ADI  Method 

The  method  proposed  by  Brian  [20]  is  similar  to  the  conventional 
ADI  method.  However,  the  successive  approximations  of  the  temperature 
are  calculated  at  the  half-time  step.  The  basic  equations  of  Brian’s 
ADI  method  are  given  as  follow: 


Ar/2  X  i,j,k  y  i,j,k  z  i,j,k 


57/2  X  i,j,k  y  i,j,k  z  i,j,k 

-V. 

=  5^u.  .  ,  +  sh.  .  ,  +  , 

At/2^  X  i,j,k  y  i,j,k  z  i,j,k 


Subtracting  Eq.  (10)  from  Eq.  (11),  we  have: 


=  sh.  .  ,  -  8^^  .  , 

S772  y  1,1, k  y  1,1, k 


Subtracting  Eq.  (11)  from  Eq.  (12),  we  have: 


^■*■1  +  ^  .  -  2V.  . 

,  -  8^^  ■  , 

5772  z  i,j,k  z  1,1, k 


(10) 

(11) 

(12) 

(11a) 

(12a) 


Eqs.  (10),  (11a)  and  (12a)  are  the  simplified  equations  suggested 
by  Brian.  After  rearranging  these  equations,  the  following  equations 
can  be  obtained: 
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! 


1 


‘^o.k  *  '^.i-l.k  ^  (13) 


■  ''i.j-l.i"  ''i.j.k  ■  ''i.M.l 

^.j-l,k-  ^.M,k 

-  t,  ,  *  +  2]  V  -  1,  < 

ijjjk'l  '■  St  J  ijjjk  i,j,k+l 

=  [^(^^.)^]  V.  .  .  +  + 

^  At  J  i,j,k  *■  At  J  iO,k  i,J»k-l  i,J,k+l 


(14) 


(15) 


Brian  shoved  that  his  scheme  is  unconditionally  stable.  However, 
there  also  exist  negative  coefficients  on  the  right-hand  sides  of  the 
discretization  equations  (13)- (15).  As  we  mentioned  earlier,  these 
negative  coefficients  are  physically  unrealistic. 


Douglas  ADI  method 

Another  unconditionally  stable  three-dimensional  ADI  method  was 
presented  by  Douglas  [19] .  The  algorithm  is  given  by  the  following 
three  equations: 


=  J  4  P’i.j.k  *  ^,j.k]  ^  ^y«io.k  *  ^z«?,j,k  (»«) 


V  - 


+  s'^d^  .  . 

z  ij,k 


(17) 
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1,1, k  1, 


^  =  i  [U-  u  +  <^  •  J  +  i  [V-  •  u  +  <^  •  ul 

2  X  L  i,j,k  ijjjk-l  2  y  >-  i,j,k  i^jk-" 
2  zL  i,j,k  i,j,kJ 


Subtracting  Eq.  (16)  from  Eq.  (17),  we  have: 


V  -  II 

i,j,k  i,.i,k  1  f2 


(17a) 


Subtracting  Eq.  (17)  from  Eq.  (18),  we  have: 


-  V 


(18a) 


Equations  (16),  (17a)  and  (18a)  are  the  simplified  equations  and 
they  can  be  rearranged  as: 


=  (^-  5)  .  ^*0.5 

*  ^,j-l,k  *  ^,j+l,k  *  ^,j,k-l  *  ^,j,k*l 


-  »■*  ’l.i-l.kMi^-l]  0.5 


Oi  j,k*  «?,j,k-  O-S  »?,j-l,k-  O-S  »?.j,l,k 


-  “•''^ik-i*  <^ik  -  o-^  «?:i.k.i 

’u.k  *  ^,j,k  -  »•=  ^,i,k.l- 
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The  unconditional  stability  of  this  algorithm  was  proved  by  Douglas 
[19] .  However,  as  in  the  Brian  ADI  method,  the  Douglas  ADI  method  has 
negative  coefficients  on  the  right-hand  sides  of  the  discretization 
equations  (19)- (21). 

New  ADI  lethod 

As  we  have  seen  above,  the  three  existing  ADI  methods  all  have 
shortcomings.  The  conventional  ADI  method  is  conditionally  stable  and 
very  small  time  steps  are  required  to  satisfy  the  stability  criterion. 
All  three  ADI  methods  have  a  common  problem:  negative  coefficients  in 
their  discretization  equations  which  are  physically  unrealistic. 

In  light  of  the  above  observation,  an  improved  ADI  method  is 
proposed.  The  conventional  three-dimensional  ADI  method  is  modified  by 
introducing  an  f- factor  (0<f<l).  Consider  a  control  volume  as  shown  in 
Fig  2.1,  the  heat  fluxes  from  the  directions  in  which  the  equation  is 
implicit  are  multiplied  by  a  factor  (3-2f)  and  the  heat  fluxes  from  the 
remaining  four  directions  are  multiplied  by  a  factor  f .  As  we  can  see, 
the  total  heat  flux  counted  in  each  direction  through  a  full  time  step 
remains  unchanged.  The  finite- difference  Eqs.  (3)- (5)  of  the 
conventional  ADI  method  are  modified  by  a  f- factor  and  become  the 
following  equations: 


=  (3-«)  *  f  ^ 


r, 


+ 


+ 


(22) 


^■*■1  -  V 

*  (»-«)  (2^) 


After  rearranging  Eq.  (22),  the  following  discretization  equation 
can  be  obtained: 


-  (^2f)  *  2(J-2f)]  .  ^  -  (3-2f) 

=  [^-  4f] 


(25) 


Similar  equations  can  be  easily  derived  from  Eqs.  (23)  and  (24)  for 
the  y-  and  z-  directions.  On  the  right-hand  side  of  Eq.  (25),  only  the 
coefficient  for  central  nodal  temperatures  at  previous  time  step  could 
be  negative.  To  avoid  it  from  becoming  negative,  we  require: 


(26) 


The  stability  criterion  can  be  determined  by  the  Von  Neumann 
method.  Assuming  that  there  exists  an  error  function  E  _  at  each 
nodal  point  in  the  following  form  [57] : 


Ep,q  r,n=  exp(i^jpAX)exp(i;52‘l^Y)exp(i)33rAZ)^“  (27) 

where  the  parameter  ^  is  the  amplification  factor  and  n  =  r/Ar,  the 
error  will  be  bounded  provided  that: 
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This  is  the  condition  for  the  solution  to  be  stable.  It  can  be 

shown  for  these  linear  problems  with  constant  coefficients  that  the 

error  function  E  ^  also  satisfies  the  finite- difference  equation 
p,q,r,n  ^ 

(25)  and  two  similar  equations  for  y-  and  z-  directions.  With  AX=AY=AZ, 
substitution  of  E  from  Eq.  (27)  into  these  equations  gives; 


3/1  -4(f)sin^(/?2AX/2)  -4(f)sin^(/?3AX/2) 
3/1  +4(g)sin^(^jAX/2) 

3/1  -4(f)sin^(;9^AX/2)  -4(f)sin^(^3AX/2) 
3/1  +4(g)sin^(^2^X/2) 


(28) 


(29) 


3/1  -4(f)sin^(^jAX/2)  -4(f)sin2(/?2AX/2) 
3/1  +4(g)sin^(^3AX/2) 


(30) 


where  g=3-2f. 

These  (2  ^3  amplification  factors  for  the 

finite- difference  equations  for  the  x- ,  y-  and  2-  directions, 
respectively.  Since  these  equations  are  used  alternately,  the  stability 
condition  should  be: 


I 

I 


^1^2^3' 


<  1 


Rearranging  (^(2^3  follows:, 


3/A  -4(f)sin^(^2^X/2)  -4(f)sin2(^3AX/2) 

[  '  '  n  ] 

^  ^  ^  3/A  +4(g)8in2(/?2AX/2) 

3/A  -4(f)8in^(/?jAX/2)  -4(f)siii^(/?3AX/2) 

3/A  +4(g)sin2(/?3AX/2)  ^ 

3/A  -4(f)8in^(^jAX/2)  -4(f)sin^(/92AX/2) 
"[ - — : — .  .  ■  - ] 


3/A  +4(g)8in^(^.AX/2) 


=  [a]«[b]«[c] 


the  stability  condition  can  be  written  as: 


t 

I 


a 


bjjcl 


<  1 


The  stability  criterion  can  be  obtained  fro«  either  one  of  the 
following  three  conditions: 

|a!<l,  |b!<l  or  |cj<l 

For  the  condition  |a|<l,  since  the  value  of  a  is  always  less  than 
unity,  we  need  only  to  consider  the  condition  a>-l.  This  leads  to: 

<  1.5 _ _ 

■  (f)8lll*(/»3AI/2)-  (g-f)sill^(^24I/2) 

It  should  be  aentioned  here  that  the  parameter  A  defined  in  Eq.  (8) 

is  always  positive.  The  right-hand  side  of  the  above  equation  has  a 

2  2 

niniauffl  value  when  sin  (^3AX/2)=1  and  sin  {^2^112)=^.  So,  the  stability 
criterion  becomes: 
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(31) 


Comparing  Eqs.  (26)  and  (31)  with  Eqs.  (7)  and  (9),  the  tine- step 
limit  for  the  conventional  ADI  method  can  now  be  increased  by  a  factor 
1/f  by  using  this  new  ADI  method.  The  computational  results  which  will 
be  discussed  later  show  that  this  modification  allows  the  time- step 
limit  to  be  increased  by  2  orders  of  magnitude  with  f=0.01  and  the 
solutions  still  remain  stable  with  high  accuracy. 

Also,  it  should  be  mentioned,  this  new  ADI  method  only  requires 
two- thirds  of  the  computer  storage  compared  to  the  Brian  or  Douglas 
method.  This  is  because  only  the  temperatures  at  one  intermediate  time 
step  need  to  be  stored. 


2.3  Results  and  Discussions 

To  validate  the  present  new  ADI  method,  the  finite- difference 
solutions  obtained  are  tested  for  a  simple  geometry  with  two  different 
boundary  conditions:  a  constant  surface  heat  flux  and  a  sudden  heating 
of  the  surface  to  a  constant  temperature.  In  addition,  it  is  compared 
with  the  Brian  and  Douglas  methods. 

Consider  a  parallelepiped  (-Lj  <  x  <  Lj,  -Lg  <  y  <  Lj,  -Lg  <  z  < 
Lg),  shown  in  Fig  2.2,  having  constant  thermophysical  properties  and 
initially  at  a  uniform  temperature  <?q=1.0.  At  time  r>0,  the 
parallelepiped  is  allowed  to  have  heat  flow  through  its  boundaries.  To 
obtain  the  temperature  distribution  within  the  parallelepiped,  equation 
(2)  must  be  solved  with  the  following  initial  conditions: 
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-1  <  X  <  1 

at  r  =  0,  <?  =  1  for  -L’  <  Y  <  (32) 

-4  <  z  <  4 


where  Lj  is  chosen  as  the  characteristic  length  4, 

L3/L1 . 


1^2—  ^  and 


Because  of  symmetry,  only  the  region  0<X<1,  O^YCL^,  and  0<Z<L2  need 
to  be  solved.  The  boundary  conditions  are: 


r  >  0 


or 


36.  36.  36.  . 

^’x-o  ^  ^\-o 

36.  _  36.  _  36. 

^  X=1  ^  Y=L’  ^  Z=L’ 

^h=l  =  ^W=L’  =  ^)z=L’  K 


i  —  'ly 

where  =  dimensionless  surface  heat  flux. 


(33a) 

(33b) 

(33c) 


In  this  report,  each  numerical  method  will  be  used  to  solve  the 
three-dimensional  heat  diffusion  equation  for  the  two  different  boundary 
conditions.  To  evaluate  the  accuracy  of  the  various  methods,  an  average 
temperature  erior  is  used.  It  is  defined  as  the  square  root  of  the 
average  of  the  squares  of  the  error  between  the  predicted  temperature 
and  the  analytical  temperature.  It  is  given  by; 


e 


I  J  K 
E  E  E 
i=l  j=l  k=l 


UK 


(34) 
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where  0^  is  the  analytical  dimensionless  temperature. 

Case  I  -  Constant  Surface  Heat  Flux 

Consider  a  parallelepiped,  initially  at  a  uniform  temperature 
<?Q=1.0.  At  time  t=0,  all  faces  of  the  parallelepiped  are  exposed  to  a 
constant  surface  heat  flux  q^=  0.5.  For  a  parallelepiped  exposed  to  a 
constant  surface  heat  flux,  the  temperature  distribution  as  a  function 
of  time  can  be  represented  by  the  summation  of  three  one- dimensional 
solutions  [17]; 


tf(X,Y,Z,T)-  Oq 


+  2q  [i 

m=0 


m=0  2^/T  2^JT 

a>  (2m+l)L’+Y  (2m+l)L’-Y 

+  E  [ierfc( - )+ierfc( - )] 

m=0  2v/t  2v/r 

0)  (2m+l)L’+Z  (2m+l)U-Z 

+  E  [ierfc( - )+ierfc( - )]} 

m=0  2v/t  2v/r 


Presented  in  Fig  2.3  are  the  results  obtained  for  a  cube  exposed  to 
a  constant  surface  heat  flux  q^=  0.5  at  time  t  -  2.0.  Twenty  nodal 
points  are  used  in  each  direction  for  this  calculation.  According  to 
Eq.  (7),  the  time- step  limit  required  for  the  conventional  ADI  method 
(f=1.0)  is  0.001875.  In  Fig  2.3,  the  solutions  from  the  conventional 
ADI  method  have  very  good  accuracy  with  time  step  0.002,  but  they  become 
unstable  as  time  step  is  increased  further.  The  Brian  and  Douglas 
methods  are  unconditionally  stable,  but  the  negative  coefficients  in  the 
discretization  equations  cause  their  solutions  to  be  physically 
unrealistic.  The  resuHs  show  that  their  solutions  have  good  accuracy 
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higiirc  2.3  Average  temperature  error  for  a  cube  with  constant 


if  time  step  is  smaller  than  0.5  but  become  more  and  more  inaccurate  if 
the  time  step  is  increased  further.  On  the  contrary,  the  proposed  new 
ADI  method  with  f  =  0.1  and  f  =  0.01  is  very  accurate  even  when  a  time 
step  of  2.0  is  used.  The  average  temperature  error  is  less  than  0.007. 
It  can  be  seen  from  Eqs.  (26)  and  (31)  that  this  f- factor  ADI  method  has 
a  much  higher  time- step  limit  than  the  conventional  ADI  method. 

Shown  in  Fig  2.4  are  the  results  at  time  r  =  10.0  for  a  cube  with 
the  same  boundary  condition  of  a  constant  surface  heat  flux  q^=  0.5. 

For  very  small  time  steps,  every  method  yields  poor  accuracy.  This  is 
due  to  the  large  amount  of  calculations  involved  and  the  accumulation  of 
round-off  errors.  For  time  steps  greater  than  0.01,  the  Brian  and 
Douglas  methods  are  always  stable  but  they  yield  poor  accuracy  with 
average  temperature  errors  up  to  about  0.15.  The  new  ADI  method  with  f 
=0.01  predicts  the  results  exceptionally  well,  the  average  temperature 
errors  are  always  less  than  0.02  for  time  steps  larger  than  0.01. 
However,  the  new  ADI  method  with  f  =  0.1  only  predicts  well  up  to  a  time 
step  of  0.5.  This  is  because  the  new  ADI  method  with  f=0.1  has  a  lover 
time  step  limit  compared  to  that  with  f=0.01. 

Figure  2.5  shows  the  variation  of  the  average  temperature  error 
with  the  f- factor  at  r  =  10.0  for  a  cube  with  the  same  boundary 
condition  of  a  constant  surface  heat  flux  q^=  0.5.  It  can  be  seen,  as 
long  as  the  solutions  do  not  diverge,  the  temperature  errors  remain 
almost  the  same  with  different  values  of  the  f- factor.  In  other  words, 
the  value  of  f  we  chose  does  not  influence  the  numerical  results  as  long 
as  the  solutions  remain  stable.  The  results  for  very  small  time  step 
A7-=0.001  always  have  larger  errors.  This  is  due  to  the  accumulation  of 
round-off  errors  we  have  mentioned  earlier.  Also,  we  can  see  that  the 
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Figure  2.4  Average  temperature  error  for  a  cube  with  constant 
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Figure  2.5  Variation  of  average  temperature  error 
constant  wall  heat  flux,  t=  10 


solutions  are  more  stable  with  smaller  values  of  the  f- factor  in  the 
sense  that  much  larger  At  can  be  used. 


Case  II  -  Constant  Vail  Temperature 

In  this  case,  the  parallelepiped,  initially  at  a  uniform 
temperature  <?q=1.0,  has  its  surfaces  suddenly  increased  and  maintained 
at  a  constant  temperature  0^  =  2.0.  The  analytical  temperature  can  be 
easily  obtained  by  using  the  method  of  separation  of  variables  [58] : 


0(X,Y,Z,r)=  1^* 


(D 

E 


CD 

s 


cos  cos  [ias#!]  cos  [IW] 


(36) 


where 


^  e^) 

T^(2m-l)(2il-lj(21-l) 


Sin 


122^111  sinM^ 


2 

'‘mnr 


^125^,2,  [i21^j2 

2  3 


Presented  in  Fig  2.6  are  the  results  obtained  for  a  cube  at  time  r 
=0.2.  At  this  time,  the  temperature  field  is  still  undergoing 
transient  development.  Similar  to  Case  I  with  constant  surface  heat 
flux,  the  conventional  ADI  method  becomes  unstable  if  the  time  step  is 
greater  than  0.002.  The  Brian  and  Douglas  methods  predict  the 
temperature  field  accurately  only  with  a  time  step  less  than  0.02. 

Their  methods  become  inaccurate  if  the  time  step  is  increased  beyond 
0.02.  The  new  ADI  method  with  both  f=0.1  and  f=0.01  always  yields 
better  accuracy  than  the  other  methods,  the  average  temperature  error 
increases  only  slightly  with  the  time  step  and  is  about  0.03  with  a  time 
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Figure  2.6  Average  temperature  error  for  a  cube  with  constant  wall  temperatuic 


step  of  0.1. 

Shown  in  Fig  2.7  are  the  results  for  a  cube  at  time  r  =  1.0.  At 
this  time,  the  temperature  field  has  already  reached  steady- state.  The 
Brian  and  Douglas  methods  predict  the  steady- state  temperature  field 
rather  poorly  if  the  time  step  is  greater  than  0.1.  The  average 
temperature  error  is  about  0.5  with  a  time  step  of  1.0.  On  the 
contrary,  the  new  ADI  method  predicts  the  steady- state  results  very 
well.  With  a  time  step  of  1.0,  the  new  ADI  method  yields  solutions  with 
an  average  temperature  error  about  0.024  for  f=0.1  and  about  0.016  for 
f=0.01. 


2.4 


Remark 


In  this  chapter,  an  f- factor  ADI  method  for  solving  transient 
three-dimensional  heat  diffusion  problems  is  introduced.  An  important 
characteristic  of  this  new  ADI  method  is  that  the  resulting 
finite- difference  equations  are  consistent  with  physical  considerations. 
Compared  to  the  conventional  ADI  method,  this  modification  allows  the 
time  step  to  be  increased  by  about  a  factor  of  1/f  without  compromising 
the  accuracy  of  the  numerical  solution.  Compared  with  the  conventional 
ADI  method  and  the  Brian  and  Douglas  ADI  methods,  this  new  ADI  method 
yields  higher  accuracy  and  requires  less  computer  storage. 
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Figure  2.7  Average  temperature  error  for  a  cube  with  constant  wall  temperatuie 


CHAPTER  3 
NUIERICAL  MODEL 


A  numerical  model  auid  solution  techniques  have  been  developed  to 
predict  the  transient  behavior  of  a  high- temperature  axially  grooved 
heat  pipe  with  thermal  energy  storage.  The  heat  transfer  through  the 
pipe  wall  and  wick,  including  the  liquid  in  the  grooves,  was  modeled  as 
three-dimensional  in  the  radial,  angular  and  axial  directions.  The 
liquid  and  vapor  flow  dynamics  were  modeled  using  quasi- steady, 
one- dimensional  methods.  The  heat  transfer  within  the  phase- change 
material,  which  was  encapsulated  in  a  cylindrical  container,  was  modeled 
as  two-dimensional  in  the  radial  and  axial  directions.  A  nodal  system 
used  to  develop  finite- difference  approximations  was  depicted  in  Figs 
3,1a  and  3.1b.  Finite- difference  equations  have  been  derived  for 
three-dimensional  heat  transfer  under  the  following  assumptions: 

(1)  The  heat  transferred  through  the  wick  and  working  fluid  is  by 
conduction  only,  since  liquid  flow  velocity  is  very  low  and 
the  liquid  thermal  conductivity  is  very  high. 

(2)  The  grooves  are  nearly  filled  with  liquid.  This  is  a  good 
assumption  for  high  temperature  heat  pipes  under  normal 
operation  without  burnout,  because  the  thermal  resistance  of 
liquid  metal  is  much  smaller  than  that  of  the  heat  pipe  wall. 

(3)  The  top  lands  of  the  groove  structure  in  the  evaporator 
section  are  adiabatic  because  no  evaporation  occurs. 

(4)  The  thermal  resistance  of  the  condensed  liquid  on  the  top 
lands  of  the  groove  structure  in  the  condenser  section  is  very 
small  compared  to  the  thermal  resistance  of  the  solid  wall  and 
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Adiabatic  Section 


3,1a  ''wlal  map  of  the  heat  pipe  -  side  view 


Nodal  map  of 


can  be  neglected. 

(5)  The  liquid- vapor  interface  temperature  is  equal  to  the  local 
vapor  temperature,  because  the  thermal  resistances  due  to 
evaporation  and  condensation  are  very  small. 

(6)  The  thermal  resistance  of  the  condensed  liquid  on  the  PCM 
cylinder  is  much  smaller  than  that  of  the  phase- change 
material,  the  surface  temperature  of  the  PCM  cylinder  is 
assumed  to  be  equal  to  the  vapor  temperature  at  the  same  axial 
location.  Hence,  only  a  two-dimensional  analysis  is  needed  to 
calculate  the  temperature  and  heat  transfer  within  the  PCM 
cylinder  because  of  angular  symmetry. 

A  variety  of  boundary  conditions  for  the  thermal  coupling  between 
the  heat  pipe  and  its  heat  source  and  sink  have  been  included  in  the 
numerical  model: 

Evaporator  surface 

(1)  variable  uniform  heat  flux,  and 

(2)  variable  uniform  temperature. 

Condenser  surface 

(1)  radiation, 

(2)  radiation  and  variable  uniform  heat  flux,  and 

(3)  radiation  and  partially- covered  variable  uniform  heat  flux. 

3.1  Heat  Conduction  Through  Pipe  Vail  and  Vick 

The  improved  three-dimensional  ADI  finite- difference  method  [18] 
was  used  to  model  the  heat  conduction  through  the  wall  and  wick. 
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including  the  liquid  in  the  grooves.  The  advantage  of  the  ADI  method  is 
that  only  tridiagonal  matrices  need  to  be  solved.  However,  the 
conventional  three-dimensional  ADI  method  is  conditionally  stable  and 
very  small  time  steps  are  required  to  ensure  convergence  and  stability. 
Since  small  Ar  is  needed  due  to  the  slender  geometry  of  the  heat  pipe, 
only  very  small  At  (about  0.001s)  can  be  used  with  the  conventional  ADI 
method.  Other  forms  of  the  ADI  method  include  the  well-known  Douglas 
[19]  and  Brian  ADI  methods [20] .  Douglas  and  Brian  ADI  methods  are 
unconditionally  stable,  they  possess  the  advantages  of  the  implicit 
scheme  with  no  limitation  on  the  size  of  time  step.  However,  Thibault 
[17]  pointed  out  that  these  two  unconditionally  stable  ADI  methods 
cannot  retain  good  accuracy  if  the  time  step  is  more  than  2  times  larger 
than  the  time- step  limit  required  for  the  conventional  ADI  method.  The 
conventional  ADI  method  was  modified  with  an  f- factor  (0<f<l)  as 
introduced  in  Chapter  2.  This  modification  allows  the  time  step  to  be 
increased  by  about  two  orders  of  magnitude  without  compromising 
significantly  on  the  accuracy  of  the  numerical  solution.  It  also  was 
shown  that  this  new  ADI  method  yields  much  higher  accuracy  than  the 
Brian  and  Douglas  ADI  methods. 

The  conventional  Douglas  and  Brian  ADI  methods  have  a  common 
problem:  negative  coefficients  in  their  discretization  equations  which 
are  physically  unrealistic  [56].  After  the  three-dimensional 
finite- difference  equations  of  conventional  ADI  method  with  cylindrical 
coordinates  are  modified  by  an  f- factor,  they  become  the  following 
equations: 


42 


(38) 


1  -i-tltjL  =  (3-2f)  ^^-pn+l/S  ^  j  +  f 

a  at/3  '  •  T  i,j,k  0  i,j,k  z  i,j,k 

™n+2/3  ,j,n+l/3 

5  -  f  *  (3-2f)  iYc^%  * '  <z’^:5{k 

|Ji+l  ^+2/3 

5  -‘^j-ays  =  f  <?^:Kk  *  *  «9T”j/3  ,  (3.2,) 


(39) 


(40) 


In  the  above  equations  the  superscripts  n,  n+1/3,  n+2/3,  and  n+1  denote 
tines  nAt,  (n+l/3)At,  (n+2/3)At,  and  (n+l)t,  respectively.  After 
rearranging  Eqs.  (38)- (40),  the  following  discretization  equations  can 
be  obtained: 


3pcr.(Ar)^  .  i  ^  /o 

(— - (3-2f)  k..(r..^)  .  (3-2f)  ki,(rj-%]  I”;V3 
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Where  p's  and  c’s  are  the  nodal  density  and  specific  heat  based  on  local 
properties  and  k’s  are  the  thermal  conductivities  based  on  the  harmonic 
mean  of  two  relevant  nodal  conductivities. 


On  the  right-hand  sides  of  Eqs.  (41)- (43),  only  the  coefficients 
for  central  nodal  temperatures  at  previous  time  step  could  be  negative. 
It  was  shown  [18]  that  the  time- step  limit  to  avoid  the  coefficients 


from  becoming  negative  can  now  be  increased  by  a  factor  of  1/f.  It  was 
also  shown  [18]  that  the  stability  criterion  can  be  increased  by  the 
same  factor.  So,  it  is  clear  that  the  time- step  limit  for  the 
conventional  ADI  method  can  now  be  increased  by  a  factor  of  1/f  by  using 
this  new  ADI  method.  The  computational  results  showed  that  this 
modification  allows  the  time- step  limit  to  be  increased  by  2  orders  of 
magnitude  with  f=0.01  and  the  solutions  still  remain  stable  with  very 
high  accuracy. 

The  resulting  finite- difference  approximations  for  all  nodes  form  a 
nonhomogeneous  set  of  linear  algebraic  equations  for  each  one- third 
time- step: 

For  the  first  one- third  time- step,: 
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for  the  second  one- third  time- step,: 
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and  for  the  third  one- third  time- step: 


Since  each  coefficient  matrix  in  Equations  (44-46)  is  tridiagonal, 
the  sets  of  equations  can  be  easily  solved  by  using  the  tridiagonal 
matrix  algorithm  (TDMA)  with  known  initial  and  boundary  conditions. 
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This  algorithm  is  very  efficient  for  digital  computation. 


3.2  Melting  and  Solidification  of  PCM 

Since  the  thermal  resistance  of  the  condensed  liquid  on  the  TES  is 
much  smaller  than  that  of  the  phase- change  material,  the  surface 
temperature  of  the  TES  is  assumed  to  be  equal  to  the  vapor  temperature 
at  the  same  axial  location.  Hence,  if  the  PCM  is  encapsulated  in  a 
cylindrical  container,  only  a  two-dimensional  analysis  is  needed  to 
calculate  the  temperature  and  heat  transfer  within  the  PCM  because  of 
angular  symmetry. 

In  the  numerical  solution  of  heat  conduction  problems  with  phase 
change,  the  heat  diffusion  equation  can  be  formulated  in  either  of  the 
following  two  ways: 


C(T)  ^  =  di*  [k(T)6rad(T)]  (47) 

or  ^  =  div  [k(H)grad(I(H))]  (48) 

Equation  (47)  is  the  basis  of  heat  capacity  methods,  while  Eq.  (48)  is 
the  basis  of  enthalpy  methods. 

In  heat  capacity  methods,  the  latent  heat  is  represented  by  a  peak 
of  small  but  finite  width  in  the  C(T)  curve  as  shown  in  Fig  3.2.  Since 
C(T)  is  not  a  continuous  function  of  Temperature,  if  a  large  time  step 
is  used  in  the  computation,  a  nodal  temperature  may  jump  past  the 
melting/freezing  temperature  range  in  one  time  step,  resulting  in  the 
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Specific  Heat,  C(T) 


latent  heat  being  ignored.  This  is  termed  "jumping  of  the  latent  heat 
peak"  and  can  be  a  major  problem.  To  avoid  it,  very  small  time  steps 
have  to  be  used. 

Enthalpy  methods  do  not  suffer  from  the  problem  of  jumping  the 
latent  heat  peak  mentioned  above  because  the  enthalpy  is  a  continuous 
function  of  temperature  as  shown  in  Fig  3.3.  However,  the  major 
disadvantage  of  enthalpy  methods  is  that  because  a  high  non  linear 
function,  T(H),  is  involved,  an  explicit  scheme  must  usually  be 
employed,  with  consequent  stability  problems.  Implicit  schemes  have 
been  proposed  by  Longworth  [59]  and  Furzeland  [60] ,  but  they  require 
iteration  at  each  time  step  and  are  less  efficient  in  terms  of  computer 
time  [33,61]. 

To  overcome  the  problem  of  jumping  the  latent  heat  peak  which  heat 
capacity  methods  suffer,  Hsiao  [31,32]  suggested  that  C(T)  should  be 
linearly  interpolated  between  the  temperatures  of  adjacent  nodes.  Hsiao 
considered  a  typical  situation  during  a  phase  change  process  (see  Fig 
3.4).  In  this  case  only  three  nodes  (inside  the  shades  region)  have 
temperatures  within  the  range  of  T^^^-AT  and  Tj^^+AT,  and  are  able  to 
include  the  latent  heat  effect  if  heat  capacity  method  is  applied.  All 
other  nodes  next  to  the  fusion  front  will  employ  the  specific  heat  of 
either  the  solid  or  liquid  phase,  depending  on  whether  the  nodal 
temperature  is  less  than  T^^j-AT  or  greater  than  T^+AT,  respectively. 
Clearly,  the  latent  heat  effect  is  not  properly  included,  and  as  a 
result  numerical  error  always  arise.  In  light  of  the  above  observation, 
Hsiao  proposed  a  new  scheme  to  improve  the  heat  capacity  model.  The 
nodal  temperature  is  not  used  directly  to  yield  the  corresponding 
specific  heat  of  the  node.  Instead,  the  following  averaged  specific 
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3.4  A  typical  location 


heat,  which  is  derived  from  the  node  (i,k)  auid  its  surrounding  nodes,  is 
employed: 


C(Ti_k)  =  ^  [C(Ti  „  , 

*  *  C(Tj  ^  ,  Tj 


(49) 


where  C(T^  ^  j^),  for  example,  represents  an  adjusted  specific 

heat,  which  is  determined  by  the  physical  status  of  the  material  with 
temperature  within  the  range  of  and  Tj,  j 

In  Hsiao’s  new  heat  capacity  method,  the  latent  heat  effect  is 
accounted  for  by  using  a  linear  interpolation  between  the  temperatures 
of  adjacent  nodes.  As  shown  in  Fig  3.4,  the  latent  heat  effect  is 
included  in  C(T^  since  the  melting/freezing  temperature 

interval  falls  in  between  ^  and  Hsiao  concluded  that  this 

new  method  can  avoid  the  problem  of  jumping  the  latent  heat  peak  and 
allows  the  use  of  a  relatively  large  time  step.  However,  we  have  tested 
Hsiao’s  new  method  but  a  large  energy  balance  error  was  found.  Pham 
[35]  also  pointed  out  that  Hsiao’s  new  method  yielded  results  with  up  to 
227,  error  compared  with  the  analytical  result.  The  low  accuracy  of 
Hsiao’s  new  method  could  be  due  to  its  ambiguous  theoretical  basis. 

Pham  [33-35]  presented  a  simple  and  accurate  method  which  combines 
the  good  features  of  enthalpy  methods  and  heat  capacity  methods.  Pham’s 
method  can  be  used  in  conjunction  with  the  two-dimensional  ADI  scheme 
using  the  following  procedure: 


1. 


♦ 

At  the  start  of  each  time  step,  the  enthalpy  change  AH  at  each 

node  is  estimated  from  the  known  temperature  ,  of  that  node  and 

1 
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those  of  its  immediate  neighbors  at  previous  time  step. 

2.  Since  the  enthalpy  is  a  continuous  function  of  the  temperature  for 


the  phase  change  material,  we  can  determine  an  estimated  new 

temperature  T.  ,  from  the  following  equation: 

1 

(50) 


where  ,  is  the  nodal  temperature  at  previous  time  step. 

1  ,K 

fq.  and  fjj  are  the  temperature  and  enthalpy  functions, 
repectively. 


3.  When  the  estimated  new  temperature  T-  ,  is  known,  the  equivalent 

1  ,K 

specific  heat  of  each  node  can  be  obtained; 


* 

c..  ,= 


=  ah 

1  ,k 

^i,k'  i,k 


(51) 


4.  With  the  equivalent  specific  heat  c-  .  known,  we  can  use  the  two 

1 

dimensional  ADI  method  to  find  the  new  nodal  temperature 

l  ,K 


One  of  the  good  features  of  Pham’s  method  is  that  it  estimates  the 

new  temperature  from  the  estimated  enthalpy  change  to  avoid  the  problem 

of  jumping  the  latent  heat  peak.  The  other  good  feature  of  Pham’s 

method  is  that  its  theoretical  basis  is  clear.  However,  Pham’s  method 

has  a  singularity  problem  in  finding  the  equivalent  specific  heat  in 

step  2.  If  there  is  no  enthalpy  change  in  a  particular  node,  the 

♦ 

estimated  new  temperature  T-  .  will  be  equal  to  the  previous  temperature 
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T?  v  Then,  we  are  not  able  to  find  the  equivalent  specific  heat  from 

1  )K 

Eq.  (51)  due  to  the  singularity  problem.  Fortunately,  we  have  found  a 
way  around  this  problem.  If  the  melting  temperature  is  and  the 
latent  heat  effect  is  over  a  2AT  interval,  let  Hj^=fjj(Tjjj- AT  ), 
H2=fH(Tjjj+AT  )  and  AH  .  We  redefine  the  equivalent  specific 

heat  in  Eq.  (13)  as  follows: 


(H2-Hj)/(2AT) 
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where  c„  and  Ci 
s  1 


are  the  specific  heats  for  solid  state  and  liquid  state. 


After  the  modification,  Eq.  (51)  is  now  used  when  only  one  of  H?  , 

1  ,K 

or  H^^j^  falls  in  between  H^  and  Hj.  In  other  words,  Eq.  (51)  can  only 
*  * 

be  used  when  AH  is  not  equal  to  zero.  Compared  with  analytical 
solutions,  this  modified  method  for  melting  and  solidification  was  found 
to  have  very  good  accuracy  and  does  not  have  the  singularity  problem 
Pham’s  method  has. 

Consider  a  one- dimensional  melting  problem  with  a  solid  in  a 
half- space  initially  at  the  melting  temperature  T^=T^=950  K.  At  time 
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zero,  the  temperature  at  the  front  surface  of  the  solid  is  suddenly 


heated  to  a  constant  Temperature  Tq=1300  K  and  melting  takes  place 
immediately.  Figs  3.5  and  3.6  show  the  comparisons  of  temperature  and 
melting  front  between  the  exact  solution  and  the  present 
finite- difference  solution  during  the  melting  process.  The  physical 
variables  chosen  in  the  numerical  experiment  are  kg=4.2  V/m-K,  kj^=2.1 
W/m-K;  Cg=6280  J/K-kg,  Cj=7370  J/K-kg;  PCM  latent  heat  is  2.58*10^ 

The  finite- difference  solution  is  obtained  by  using  100  equally  spaced 
elements  along  a  length  1=0.01  m.  Time  step  At=0.01s  and  a  phase 
transition  temperature  interval  2AT=  2.0  K  are  used  for  this  example. 

The  dimensionless  time  is  defined  as  r=t/(l  /oj) .  The  Stefan  number  is 
defined  as  Ste=C2(TQ-T^)/L  and  is  equal  to  unity  in  this  problem.  As  it 
can  be  seen  in  Figs  3.5  and  3.6,  the  present  solution  agrees  very  well 
with  the  exact  solution. 


3.3  One- Dimensional  Vapor  Flow  Model 

The  vapor  flow  was  modeled  by  using  a  quasi- steady ,  one- dimensional 

friction  coefficient  developed  by  Bowman  [46] .  In  the  evaporation 

region,  mass  blowing  causes  a  slight  steepening  in  the  velocity 

gradients  at  the  pipe  wall,  leading  to  an  increase  in  the  friction 

coefficient.  Bowman  pointed  out  that  the  favorable  pressure  gradient  in 

the  mass  blowing  region  influences  the  flow  to  remain  laminar,  even  for 

fi 

very  large  axial  Reynolds  numbers  up  to  10°.  In  the  condenser  region, 
where  there  is  mass  removal  and  an  adverse  pressure  gradient,  the  flow 
was  found  to  stay  laminar  at  axial  Reynolds  numbers  around  12000.  In 
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this  research,  the  vapor  flow  was  always  assumed  to  be  laminar  because 
the  maximum  axial  Reynolds  numbers  are  always  much  less  than  12000.  The 
correlation  of  the  vapor  friction  coefficient  for  laminar  flow  given  by 
Bowman  can  be  expressed  as: 

f  =  (1.2337-0.2337  ’^^w)  (53) 


where  Ma  is  the  Mach  number  based  on  the  local  mean  axial  velocity  U, 
Re^  is  the  radial  Reynolds  number,  and  Re  is  the  axial  Reynolds  number 
defined  as: 


Re= 


In  these  expression,  p  is  the  vapor  flow  density,  /i  is  vapor  dynamic 
viscosity,  v  is  the  radial  velocity  at  the  wall  and  is  the  hydraulic 
diameter  of  the  vapor  core. 


The  vapor  flow  was  assumed  to  be  compressible,  one- dimensional  and 
quasi- steady.  The  governing  equations  for  such  a  flow  can  be  expressed 
in  terms  of  influence  coefficients  as  presented  by  Shapiro  [62]: 


=  Fc  4f 


dz  dib 

”i;  ^ 


(54) 


with  the  two  influence  coefficients  given  by: 


f  ,a' 


4Ma^[l+  ^  Ma^] 


1-  Ma" 
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(56) 


2[l+7lla^]  [1+  ^  Ha^] 


A, a 


1-Ma" 


where  f  is  the  friction  coefficient  defined  earlier  by  Eq.  (53),  7  is 
the  ratio  of  specific  heats,  z  is  the  axial  coordinate,  and  A  is  the 
mass  flow  rate. 

For  the  friction  solution,  a  second  expression  is  needed  to  relate 
the  change  in  total  pressure  (P^)  to  the  change  in  mass  flow  rate  and  to 
the  friction  coefficient.  From  Shapiro  [62]: 


dP 


®  =  Fx  .4f  +  F. 


dih 
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Other  useful  relations  are: 
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These  relate  propei’ties  between  two  different  axial  locations,  a 
and  b,  in  the  vapor  flow  field. 

3.4  Coupling  of  Vapor  Flow  With  Evaporation  and  Condensation 

To  calculate  the  pressure  and  temperature  variations  in  the  vapor 
flow,  we  need  to  know  the  evaporation  and  condensation  rates.  However, 
these  rates  depend  on  the  vapor  temperature  and  liquid- vapor  interface 
temperature  distributions.  This  means  that  the  vapor  pressure  and 
temperature  variations  are  coupled  to  evaporation  and  condensation  rates 
simultaneously.  In  the  present  model,  the  evaporation  and  condensation 
rates  are  coupled  to  the  vapor  temperature  and  pressure  in  an  explicit 
manner  so  that  no  iterations  are  required.  However,  we  still  have  to 
guess  the  vapor  temperature  at  the  evaporator  end  to  calculate  the  vapor 
temperature  distribution.  The  vapor  temperature  at  the  evaporator  end 
can  be  estimated  based  on  an  assumption.  Since  the  vapor  density  is 
very  small,  we  can  assume  that  the  heat  absorbed  by  the  vapor  itself  is 
negligible  compared  to  the  total  evaporation  and  condensation  rates.  In 
other  words,  at  every  time  step,  the  total  evaporation  rate  and  the 
total  condensation  rate  are  equal.  So,  the  following  equation  should  be 
always  satisfied: 

where  is  the  total  evaporation  rate  and  is  the  total 
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condensation  rate  (see  Fig  3.7). 

In  this  research,  the  one- dimensional  vapor  flow  was  coupled  with 
the  local  evaporation  and  condensation  rates,  as  shown  in  Fig  3.7,  using 
the  following  procedure; 

(1)  The  local  evaporation  rate  or  condensation  rate  Qy(k)  at  previous 
time  step,  can  be  evaluated  from  previous  local  temperatures 

(2)  Guess  the  new  vapor  temperature  at  the  evaporator  end  T^^^(l),  and 
then  use  the  evaporation  and  condensation  rates  to  predict  the 
vapor  pressure  and  temperature  variations  at  new  time  step  (Py^^(k) 
and  T;^‘(k)). 

(3)  As  T;*‘(k)  is  obtained,  we  can  calculate  the  new  temperature 

distribution  in  the  heat  pipe  wall  and  wick  . 

(4)  The  local  evaporation  or  condensation  rate  ^(k)  at  new  time  step 
can  be  evaluated  from  T^^^(k)  and 

(5)  The  total  evaporation  rate  £  and  condensation  rate  E  can 
be  evaluated  from  Qy^^(k). 

(6)  If  A=  E  -  E  >  t,  go  back  to  step  (2)  and  iterate  until  A 
<  c. 

(7)  After  it  converges,  Qy^^(k),  T^'^^(k),  Py''^(k),  and  T^^j^  at  new 
time  step  can  finally  be  obtained. 
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Figure  3.7  Coupling  of  vapor  flow  with  evaporation  and  condensation 
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3.5  Liuid  Pressure  Drop 

In  calculating  the  liquid  pressure  drop,  we  assumed  that  the  groove 
structures  are  fully  wetted  with  working  fluid  and  the  liquid  flow  is 
always  laminar  owing  to  the  generally  low  liquid  velocity.  The  liquid 
pressure  drop  can  be  obtained  by  using  the  following  expression  [63] : 

dPi 

32"  -  -FlQ  (64) 

Here,  Q  is  the  local  axial  heat  rate  along  a  groove,  and  is  a 
frictional  coefficient  for  the  liquid  flow  and  is  defined  as: 

h  =  Egbr  <“> 

6  1 

Here,  pj  is  the  liquid  viscosity,  is  the  groove  cross-sectional  area, 

X  is  the  latent  heat  of  evaporation,  p-^  is  the  liquid  density,  and  K  is 
the  permeability  of  the  groove  structure  and  is  calculated  from  the 
equation: 

2cr?  , 

*  =  (IjHep 

where  c  is  the  groove  porosity,  r^^  ^  is  the  hydraulic  radius  of  the 
groove  structure  defined  as  twice  the  cross-sectional  area  divided  by 
the  wetted  perimeter,  and  (f^Re^)  is  a  constant  for  laminar  flow  whose 
magnitude  depends  only  on  the  geometry  of  the  groove  structure  and  can 
be  obtained  from  Ref.  63. 


63 


CHAPTER  4 


TRANSIENT  BEHAVIOR  OF  HP/TES  SYSTEM 
TINDER  PULSE  HEAT  LOADS  APPLIED  AT  THE  EVAPORATOR 


In  this  chapter,  we  will  examine  two  important  points.  First,  how 
effectively  can  a  PCM  mitigate  the  adverse  effects  of  pulse  heat  loads 
applied  at  the  evaporator?  Second,  is  there  any  disadvantage  associated 
with  installing  phase- change  material  (PCM)  in  the  vapor  core?  Also,  if 
there  are  significant  difficulties  associated  with  such  a  system,  will 
they  be  offset  by  the  capabilities  of  the  PCM  itself? 

It  is  clear  that  the  main  disadvantage  of  installing  phase  change 
material  in  the  vapor  core  of  a  heat  pipe  is  the  accompanying  reduction 
in  the  vapor  flow  area.  This  reduction  in  vapor  flow  area  could  cause 
vapor  pressure  drop  and  vapor  velocity  to  increase,  and  thus  decrease 
the  heat  pipe  capability.  Fortunately,  the  PCM  itself  can  absorb  a 
large  portion  of  the  heat  loads  during  the  melting  process  after  pulse 
heat  loads  are  applied.  Some  vapor  will  condense  on  the  surface  of  PCM 
cylinders  and  reduce  the  vapor  velocity  and  pressure  drop.  The  net 
increase  (or  decrease)  in  vapor  velocity  and  pressure  drop  due  to  the 
installation  of  PCM  is  strongly  dependent  on  how  efficiently  the  PCM  can 
absorb  the  heat  loads. 


4.1  Limitations  of  the  Heat  Pine 


Circulation  of  working  fluid  is  an  important  heat  pipe  design 
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factor.  The  greatest  possible  circulation  is  required  to  obtain  a 
maximum  heat  transport  capability  from  the  heat  pipe.  Limitations  on 
the  heat  transport  capability  include  capillary  pumping  ability 
(capillary  limit),  choking  of  vapor  flow  (sonic  limit),  tearing  of 
liquid  off  the  liquid- vapor  interface  by  vapor  flowing  at  high  velocity 
(entrainment  limit),  and  disruption  of  the  liquid  flow  by  nucleate 
boiling  in  the  wick  (boiling  limit). 

Figure  4.1  shows  how  transport  capability  varies  with  the  operating 
temperature  for  a  screen- wrapped  grooved  heat  pipe  without  anything 
installed  in  the  vapor  core.  These  data  are  based  on  a  total  heat  pipe 
length  of  1.0  m  with  evaporator,  condenser  and  adiabatic  sections  of 
0.3ra,  0.3m  and  0.4m,  respectively.  The  heat  pipe  outside  diameter  is 
assumed  to  be  1.9  cm  (3/4  in)  with  an  inside  diameter  of  1.4  cm.  Liquid 
sodium  is  used  as  the  working  fluid,  and  the  screen  wick  has  200  meshes 
per  inch.  As  one  can  see,  within  the  operating  temperature  range 
between  900  and  1300  K,  the  entrainment  limit  places  the  greatest 
restrictions  on  heat  transport  capability.  At  950  K,  the  entrainment 
limit  is  about  5.4  kV,  which  is  equivalent  to  a  uniform  heat  flux  30 
W/cm  applied  at  the  evaporator.  If  an  empty  cylinder  was  installed  in 
the  vapor  core,  the  heat  pipe  capacity  would  be  degraded  due  to  the 
reduction  in  vapor  flow  area.  Fig  4.2  shows  the  operation  limits  for 
the  heat  pipe  with  an  empty  cylinder  with  a  radius  of  0.4  cm  mounted  in 
the  vapor  core.  The  entrainment,  capillary,  and  sonic  limits  are 
decreased  due  to  installation  of  the  empty  cylinder.  At  an  operating 
temperature  of  950  K,  the  entrainment  limit  is  decreased  by  about  337i  to 
3.6  kV,  which  is  equivalent  to  a  uniform  heat  flux  20  V/cm^  applied  at 
the  heat  pipe  evaporator.  Fortunately,  filling  the  cylinder  with  PCM 
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Figure  4.1  Operation  limits  for  a  grooved  heat  pipe  without  anything  in.stalled 
in  the  vapor  core 
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Operation  limits  for  a  grooved  heat  pipe  with  an  empty  eylind 
mounted  in  the  vapor  core 


can  offset  this  disadvantage  by  absorbing  a  large  portion  of  the  heat 
loads  during  the  melting  process.  Such  a  mechanism  will  be  discussed  in 
the  last  section  of  this  chapter. 


4.2  Operation  Under  Normal  Conditions 

In  this  section,  we  will  see  how  the  HP/TES  system  responds  under 
the  pulse  heat  loads  which  are  less  than  the  heat  pipe  operation  limit. 
We  have  applied  the  numerical  model,  introduced  in  Chapter  3,  to  a  heat 
pipe  with  18  grooves  using  sodium  as  the  working  fluid.  The  material 
for  the  heat  pipe  container  is  assumed  to  be  type  316  stainless  steel. 
All  axial  and  radial  dimensions  of  the  heat  pipe  are  the  same  as  those 
just  discussed  in  connection  with  Fig  4.1.  A  uniform  heat  flux  is 
applied  to  the  evaporator,  and  heat  is  removed  at  the  condenser  by 
radiative  heat  transfer  to  the  surroundings  which  are  maintained  at  0  K. 

The  eraissivity  of  the  condenser  wall  surface  is  assumed  equal  to  unity. 

R 

Lithium  hydride,  which  has  a  latent  heat  of  2.58*10  J/Kg  and  a  melting 
temperature  of  956  K,  is  used  as  the  phase- change  material. 

For  numerical  modeling  of  the  heat  pipe  wall  and  wicks,  8  and  40 
nodes  were  chosen  in  radial  and  axial  directions,  respectively.  Since 
the  heat  pipe  could  be  divided  into  18  identical  land- and- groove 
subunits,  only  4  nodes  were  needed  in  the  angular  direction  for  each 
unit  once  symmetry  was  invoked.  The  transient  responses  of  three 
different  HP/TES  configurations  were  compared:  (1)  a  heat  pipe  with  one 
large  empty  cylinder  installed  in  the  vapor  core,  (2)  a  heat  pipe  with 
one  large  PCM  cylinder,  and  (3)  a  heat  pipe  with  six  small  PCM 
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cylinders.  The  radius  of  the  single  large  cylinders  is  0.4  cm  while 
each  of  the  small  cylinders  has  a  radius  of  0.163  cm.  These  radii  were 
chosen  so  that  the  large  cylinder  holds  the  same  amount  of  PCM  as  the 
six  small  cylinders  combined.  The  hydraulic  diameter  for  vapor  flow  is 
about  0.82  cm  for  the  heat  pipe  with  a  single  large  cylinder  and  about 
0.75  cm  for  the  one  with  six  small  PCM  cylinders.  For  the  numerical 
modeling  of  the  PCM,  40  nodes  in  the  radial  direction  were  chosen  for 
the  large  PCM  cylinder  and  16  nodes  were  used  with  the  small  one.  A 
phase  transition  temperature  interval  of  2^1^=  10.0  K  was  assumed  and 
the  time  step  At=  0.1  s  has  been  used  for  all  the  examples  in  this 
report . 

Figure  4.3  illustrates  the  transient  response  of  four  different 
HP/TES  configurations  when  a  higher  heat  load  is  suddenly  applied  to  the 
evaporator.  Before  t=  10  s,  all  four  different  heat  pipes  are  operating 
at  steady- state  under  a  uniform  evaporator  heat  load  q=  4.3  V/cm  .  The 
average  heat  pipe  temperature  is  about  940  K.  After  t=  10  s,  a  higher 
heat  load  of  q=  10  V/cm  is  suddenly  applied  to  the  evaporator.  As 
shown  in  the  figure,  the  temperature  of  both  heat  pipes  without  PCM 
increases  quite  sharply.  It  should  also  be  noted  that  the  temperature 
of  the  heat  pipe  without  anything  installed  in  the  vapor  core  shows  no 
difference  with  that  of  the  one  mounted  with  an  empty  cylinder  through 
the  whole  test  period.  The  temperature  of  both  heat  pipes  fitted  with 
PCM  also  increases  rapidly  immediately  right  after  the  heat  load  is 
applied,  but  this  steep  temperature  increase  is  arrested  when  the  PCM 
reaches  its  fusion  temperature  where  melting  begins.  One  should  also 
note  that  the  temperature  increase  of  the  heat  pipe  with  six  small  PCM 
cylinders  is  much  slower  than  for  the  one  with  a  single  large  PCM 
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Figure  4.3  Transient  response  of  heat  pipes  with  a  sudden  increase  in  heat  input 
from  q=  4.3  to  10  W/cm^ 


cylinder  during  the  melting  process-  The  smaller  PCM  cylinders  have  an 
advantage  because  their  total  heat  transfer  surface  area  is  greater  and 
their  heat  conduction  path  is  shorter  compared  to  the  configuration  with 
a  single  large  PCM  cylinaer.  However,  the  six  small  PCM  cylinders  will 
be  completely  melted  earlier  than  a  single  large  PCM  cylinder.  The  six 
small  PCM  cylinders  are  completely  melted  at  about  t=  100  s,  and  the 
larger  PCM  cylinder  does  not  completely  melt  until  about  t=  130  s. 

After  the  PCM  has  completely  melted,  the  temperature  of  both  heat  pipes 
with  PCM  starts  to  increase  rapidly  until  the  heat  pipes  reach  a  new 
steady-  state  condition  with  q=  10  V/cm  .  Both  heat  pipes  without  PCM 
reaches  their  new  steady- state  conditions  at  about  t=  350  s.  Due  to  the 
additional  heat  capacity  of  the  phase- change  materials,  the  other  two 
heat  pipes  with  PCM  reach  the  steady- state  condition  somewhat  later  at 
about  t=  580  s.  Ve  also  applied  a  lumped  model  to  predict  the  heat  pipe 
transient  behavior.  As  one  can  see,  the  results  from  the  lumped  model 
for  the  heat  pipe  without  a  PCM  are  in  good  agreement  with  those  from 
the  finite- difference  method.  The  heat  pipe  temperature  predicted  by 
the  lumped  model  averages  about  10  K  average  lower  than  that  obtained 
using  the  finite- difference  method.  This  discrepancy  arises  because  the 
heat  removed  from  the  condenser  by  radiation  is  overestimated  by  the 
lumped  model,  which  uses  the  average  heat  pipe  temperature  as  the 
condenser  wall  surface  temperature. 

An  interesting  output  from  the  finite- difference  solution  is  the 
fraction  of  the  heat  conducted  through  the  heat  pipe  wall  at  the 
evaporator  which  is  absorbed  by  the  PCM  during  the  melting  process.  In 
other  words,  one  would  like  to  know  what  percentage  of  the  heat 
conducted  into  the  vapor  core  significantly  contributes  to  the  vapor 
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flow.  As  shown  in  Fig  4.4,  we  define  Q-  as  the  total  heat  rate 

transferred  through  the  heat  pipe  wall  at  the  evaporator  and  Q.  as  the 

total  heat  rate  absorbed  by  the  PCM.  Both  0-  and  Q.  are  estimated  from 
the  temperature  distributions  near  the  heat  pipe  inner  wall  and  the 
outside  surfaces  of  the  PCM  cylinders. 

Figure  4.5  plots  the  ratio  between  Q.  and  Q-  for  the  same 

u  lU 

operating  conditions  depicted  in  Fig  4.3.  For  the  heat  pipe  with  a 
single  large  PCM  cylinder,  this  ratio  reaches  467,  and  then  declines  to 
277i  before  the  PCM  is  fully  melted.  If  the  heat  pipes  is  equipped  with 

six  small  PCM  cylinders,  the  ratio  of  to  reaches  527.  and  then 
decreases  to  407.  before  the  PCM  is  fully  melted.  Apparently  if  the  heat 
pipe  contains  a  single  large  PCM  cylinder,  no  more  than  737.  maximum  of 
the  heat  conducted  through  the  heat  pipe  wall  co'  tj"’ites  to  the  vapor 
flow  during  the  PCM  melting  process.  For  the  TES  configuration  using 
six  small  PCM  cylinders,  the  maximum  fraction  influencing  vapor  flow  is 
only  b07..  These  results  imply  that  the  increases  in  vapor  pressure  drop 
and  vapor  velocity  due  to  installation  of  the  PCM  can  be  significantly 
offset  by  the  capabilities  of  the  PCM  itself  during  the  melting  process. 

Figures  4.6a  and  4.6b  show  the  axial  variation  of  vapor  pressure 
and  temperature  for  four  different  HP /TES  configurations  at  t=  100  s. 

The  heat  pipes  without  PCM  have  vapor  temperature  of  1082  K,  which  is 
much  higher  than  that  of  the  other  two  heat  pipe  designs  which 
incorporate  PCM.  The  heat  pipe  with  a  single  large  cylinder  of 
phase- change  material  has  a  vapor  temperature  of  985  K,  and  distributing 
the  PCM  among  six  small  cylinders  lowers  the  vapor  temperature  to  only 
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Figure  4.6a  Axial  variation  of  vapor  pressure  at  t= 
conditions  depicted  in  Fig.  4.3 


968  K.  The  results  in  Figs  4.6a  and  4,6b  indicate  that  the  vapor 
pressure  and  temperature  drops  are  strongly  dependent  on  the  operating 
vapor  temperature.  A  lower  vapor  temperature  will  result  in  a  lower 
speed  of  sound  in  the  vapor  and  a  higher  Mach  number  for  the  flow.  From 
Eq.  (53)  in  Chapter  3,  one  can  see  that  higher  Mach  numbers  will  lead  to 
higher  friction  coefficients,  resulting  in  larger  total  pressure  drops. 
Higher  Mach  numbers  also  cause  the  static  pressure  to  drop  even  more. 

In  the  adiabatic  section,  much  less  mass  blowing  or  suction  will  take 
place.  Thus,  the  pressure  decreases  there  are  mainly  due  to  friction 
and  will  be  smaller  than  those  occurring  in  the  evaporator.  As 
expected,  these  pressure  and  temperature  losses  are  partially  recovered 
in  the  condenser  section.  Fig  4.6a  also  shows  the  vapor  pressure  drop 
of  the  heat  pipe  without  anything  installed  in  the  vapor  core  is  smaller 
than  that  of  the  one  mounted  with  an  empty  cylinder  due  to  larger  vapor 
flow  area. 

In  Fig  4.7,  the  results  from  Fig  4.3  are  compared  with  those 
obtained  using  a  larger  time  step  and  wider  grid  spacings.  The  time 
step  was  increased  from  0.1  second  to  0.5  second  and  the  grid  spacings 
Ar,  Ar^  and  Az  have  been  doubled.  As  shown  in  the  figure,  the  longer 
time  step  and  coarser  grid  mesh  have  little  effect  on  the  solution.  For 
all  practical  purposes,  the  numerical  results  presented  in  Fig  4.3  are 
essentially  independent  of  the  time  step  and  grid  spacing. 

The  transient  responses  of  the  heat  pipes  with  a  pulse  heat  load 
applied  to  the  evaporator  from  t=  10  s  to  t=  100  s  are  shown  in  Figs  4.8 
and  4.9.  As  can  be  observed  from  Fig  4.8,  the  heat  pipe  responds  very 
quickly,  and  the  temperature  starts  to  decrease  as  soon  as  the  pulse 
heat  load  is  removed  at  t=  100  s.  The  temperature  of  the  heat 
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Figure  (4.7)  Comparison  of  results  from  Fig.  4.3  and  thos^  obtained  using  a 
larger  time  step  and  wider  grid  spacings 


1100- 


Figure  4.8  Transient  response  of  heat  pipes  with  pulsed  heat  loads  (p  10  W/em 


w/  TES,  Nt=1,  Ri=0.4  cm 


Figure  4.{  Portion  of  PCM  melted  versus  time  for  pulsed  heat  loads  q=  10  W/cm 


pipe  without  a  PCM  decreases  quite  rapidly  after  the  pulse  heat  load  is 
removed.  The  temperature  of  the  other  two  designs  with  PCM  also 
decreases  immediately  after  the  pulse  heat  load  is  removed,  but  the 
decrease  becomes  very  slow  when  the  PCM  reaches  its  melting  point  and 
starts  to  solidify.  The  six  small  PCM  cylinders  will  become  completely 
solidified  sooner  than  a  single  large  PCM  cylinder.  After  the  PCM 
cylinders  have  completely  solidified,  the  temperature  of  both  heat  pipes 
incorporating  a  PCM  begins  decreasing  again  and  the  heat  pipes  return  to 
their  initial  steady- state  condition  with  q=  4.3  V/cm  .  For  the  heat 
pipe  with  six  PCM  cylinders  and  the  design  with  one  large  PCM  cylinder 
to  return  to  the  initial  steady- state  conditions  after  the  pulse  heat 
load  is  removed  requires  about  1,620  seconds  and  1,740  seconds, 
respectively.  Fig  4.9  shows  the  percentage  of  PCM  melted  versus  time 
for  the  same  cases  covered  in  Fig  4.9.  One  can  see  that  the 
phase- change  material  does  not  respond  as  fast  as  the  heat  pipe 
temperature  does.  In  fact,  after  the  pulse  heat  load  is  removed  at  t= 
100  s,  27i  more  of  the  six  small  PCM  cylinders  and  107i  more  of  the  single 
large  PCM  cylinder  will  still  be  melted  by  residual  heat.  It  should 
also  be  noted  that  the  small  PCM  cylinders  solidify  more  quickly  than 
the  larger  one  can. 

Figures  4.10  and  4.11  illustrate  the  results  for  the  transient 
response  of  the  heat  pipes  with  periodic  pulse  heat  loads.  The  time 
period  of  the  pulses  is  200  s  and  each  pulse  heat  load  lasts  for  20  s. 
These  periodic  pulse  heat  loads  are  applied  before  the  heat  pipes  have 
enough  time  to  return  to  their  initial  steady- state  operating  condition. 
The  temperature  response  on  each  cycle  is  similar  to  the  results  shown 
in  Fig  4.8.  On  the  other  hand,  temperature  of  the  heat  pipe  without 
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Figure  4.U  Transient  response  of  heat  pipes  with  periodic  pulsed  heat  load 
0=  10  V/c«2 
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Figure  4.11  Portion  of  PCM  melted  versus  time  for  periodic  pulsed  heat  loads 


a  PCM  temperature  regulation  system  goes  up  and  down  periodically.  The 
temperature  of  the  heat  pipe  with  six  small  PCM  cylinders  remains  almost 
constant  over  several  cycles  due  to  very  efficient  melting  and 
solidification.  It  is  clear  from  Fig  4.11  that  the  percentage  of  molten 
PCM  for  both  heat  pipes  containing  PCM  cylinders  continues  to  increase 
as  the  cycles  persist.  This  increase  in  the  melted  PCM  fraction  over 
time  takes  place  because  the  PCM  does  not  have  enough  time  during  the 
low- power  periods  to  solidify  completely  before  the  next  pulse  cycle 
begins. 


4.3  Operation  Near  HP /TES  Limitation 

As  was  mentioned  earlier,  the  main  drawback  to  installing  PCM  in 
the  vapor  core  of  a  heat  pipe,  neglecting  the  PCM  ability  to  absorb  a 
portion  of  the  heat  load  during  melting,  is  that  the  heat  pipe 
capability  is  degraded  due  to  a  reduction  in  the  vapor  flow  area.  From 
Fig.  4.2,  the  heat  transport  limitation  for  a  heat  pipe  with  a  large 
empty  cylinder  at  temperature  950  K  is  equivalent  to  a  uniform  heat  flux 
of  20  V/cra  applied  to  the  evaporator.  In  this  section,  we  will  examine 
the  transient  response  of  the  HP/TES  system  under  pulse  heat  loads  near 
this  heat  pipe  limitation  to  see  if  the  decrease  in  heat  pipe  capability 
from  the  reduction  in  vapor  flow  area  can  be  recovered. 

Figure  4.12  shows  the  transient  response  of  three  different  HP/TES 
configurations  when  a  pulse  heat  load  near  the  heat  pipe  limitation  is 
suddenly  applied  to  the  evaporator.  Before  t=  10  s,  all  three  different 
heat  pipes  are  operating  at  steady- state  under  a  uniform 
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Figure  4.12  Transient  response  of  heat  pipes  with  a  sudden  increase  in  heat 
inputs  from  q=  4.3  to  20  W/cm^ 


2 

heat  load  q=  4.3  V/cm  in  the  evaporator.  The  average  heat  pipe 
temperature  is  about  940  K.  After  t=  10  s,  a  higher  heat  load  of  q=  20 
V/cm  will  be  applied  to  the  evaporator.  As  shown  in  the  figure,  the 
trend  of  the  temperature  increases  for  all  three  heat  pipes  is  similar 
to  the  results  seen  in  Fig  4.3  with  lower  pulse  heat  loads.  However, 
the  PCM  will  be  fully  melted  much  earlier  under  the  higher  pulse  heat 
load.  The  six  small  PCM  cylinders  will  finish  melting  at  about  t-  50  s 
while  the  single  larger  PCM  cylinder  does  not  completely  melt  until 
about  t=  70  s. 

Figure  4.13  depicts  the  axial  variation  of  liquid  and  vapor 
pressure  for  the  heat  pipe  without  a  PCM  at  t=  12  s.  At  this  moment, 
the  heat  pipe  has  the  greatest  liquid  and  vapor  pressure  drops  due  to 
the  high  liquid  and  vapor  mass  flow  rates  and  low  heat  pipe  temperature. 
Apparently,  the  vapor  pressure  drop  dominates  the  liquid  pressure  drop 
for  this  type  of  heat  pipe.  As  one  knows,  the  heat  pipe  capillary  limit 
is  strongly  dependent  on  the  overall  liquid  and  vapor  pressure  drops. 
Thus,  for  this  type  of  heat  pipe,  the  effect  of  the  liquid  pressure  drop 
on  the  heat  pipe  capillary  limit  can  be  neglected. 

The  total  heat  rates  transferred  through  the  heat  pipe  wall  at  the 
evaporator  for  three  different  heat  pipes  are  plotted  in  Fig  4.14.  The 

total  pulse  heat  load  applied  at  evaporator  is  Qg=  3.64  kV  (which  is 
equivalent  to  a  uniform  heat  flux  q=  20  V/cm  )  after  t=  10  s.  For  the 

heat  pipe  without  a  PCM,  only  2.8  kV  is  transferred  through  the 
heat  pipe  wall  into  the  vapor  core  at  t=  20  s.  This  reduction  in  the 
heat  transfer  rate  arises  because  a  large  portion  of  the  heat  load  is 
absorbed  by  the  heat  pipe  wall  due  to  the  rapid  temperature  increases 
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Figure  4.13  Axial  variation  of  liquid  and  vapor  pressure  at  t=  12  s  for  the  same 
operating  conditions  depicted  in  Fig.  4.12 


Figure  4,14  Variation  of  Q-  for  the  same  operating  conditions  depicted  in 


taking  place  at  this  moment.  The  value  rises  gradually  once  the 

heat  pipe  temperature  increase  slows  down,  and  finally  approaches 
upon  reaching  another  steady- state.  On  the  other  hand,  for  both  heat 
pipes  with  PCM  most  of  the  heat  loads  are  transferred  through  the  heat 
pipe  wall  into  the  vapor  core  at  t=  20  s.  Heat  transfer  to  the  vapor 
core  is  more  efficient  for  the  PCM  designs  because  the  heat  pipe 
temperature  increases  are  being  slowed  by  the  PCM  melting  process.  The 

value  for  both  heat  pipes  with  PCM  remains  almost  constant 
throughout  the  melting  process.  It  then  drops  suddenly  as  soon  as  the 
PCM  is  completely  melted,  and  the  heat  pipe  temperature  starts  to 
increase  rapidly. 

Figure  4.15  shows  the  ratio  for  both  heat  pipes  equipped 

with  PCM  and  operating  under  pulse  heat  loads  near  the  heat  pipe 
limitation.  It  is  quite  apparent  that  the  PCM  responds  very  quickly  to 
pulse  loading.  The  PCM  starts  to  melt  and  absorbs  a  large  portion  of 
the  heat  immediately  after  the  pulse  heat  loads  are  applied.  For  the 
heat  pipe  with  one  large  PCM  cylinder,  this  ratio  reaches  a  maximum  of 
587i  and  then  decreases  to  387t  before  the  PCM  is  completely  melted.  For 

the  one  with  six  small  PCM  cylinders,  the  ratio  peaks  at  717.  and 

then  declines  to  537.  when  the  PCM  is  totally  melted.  In  effect,  the 
heat  pipe  with  one  large  PCM  cylinder  has  at  most  627.  of  the  heat  which 
is  conducted  through  the  heat  pipe  wall  contributing  to  the  vapor  flow 
during  the  PCM  melting  process.  In  the  design  with  six  small  PCM 
cylinders,  the  maximum  fraction  adding  to  the  vapor  flow  is  only  477.. 

As  we  mentioned  earlier,  the  decrease  in  the  vapor  flow  area  due  to  the 
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installation  of  PCM  causes  the  heat  pipe  capability  to  decrease  by  about 
337*.  From  the  results  shown  in  Fig  4.15,  it  is  clear  that  this  decrease 
in  heat  pipe  capability  can  be  completely  offset  by  the  capabilities  of 
the  PCM  itself  during  the  melting  process. 

Comparing  the  results  between  Fig  4.5  and  Fig  4.15,  one  finds  that 

the  ratio  Q+/Q-_  is  larger  with  higher  pulse  heat  loads  during  the  PCM 
melting  process.  Since  the  heat  pipe  temperature  remains  about  the  same 
during  the  PCM  melting  process  for  both  pulse  heat  loads,  the  heat 
removed  from  the  condenser  by  radiation  will  also  be  similar.  Thus,  the 
percentage  of  the  heat  loads  directly  contributing  to  the  vapor  flow 
during  the  PCM  melting  process  will  be  lower  under  higher  pulse  heat 
loads.  Ve  can,  therefore,  predict  that  the  HP/TES  system  can  operate 
without  burnout  under  uniform  heat  fluxes  even  greater  than  q=  30  V/cm  , 
the  heat  pipe  limitation  present  without  anything  installed  in  the  vapor 
core.  Not  only  can  the  PCM  recover  the  decrease  in  heat  pipe 
performance  due  to  the  reduction  in  vapor  flow  area,  but  it  can  actually 
increase  heat  pipe  transport  capability. 

In  the  design  of  a  HF/TES  system,  one  should  choose  a  PCM  with  a 
melting  point  slightly  higher  than  the  normal  operating  temperature. 

Then  if  a  pulse  heat  load  higher  than  the  heat  pipe  limitation  is 
applied,  the  PCM  can  respond  fast  enough  to  begin  melting  and  absorb 
some  of  the  heat  before  the  heat  pipe  reaches  its  operating  limit  and 
burns  out.  To  reduce  the  chance  of  completely  melting  during  the  pulse 
period,  the  latent  heat  of  fusion  of  the  chosen  PCM  should  also  be  as 
large  as  possible. 

The  concept  of  incorporating  phase- change  material  inside  a 
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low- temperature  heat  pipe  (such  as  a  water  heat  pipe)  is  also  sound  if 
the  goal  is  to  limit  the  temperature  extremes  encountered  when  the  heat 
load  is  time- dependent.  For  most  low  temperature  heat  pipes,  the  vapor 
pressure  drop  is  small,  and  vapor  flow  usually  does  not  play  an 
important  role  in  determining  the  heat  pipe  capability.  Thus  the 
increases  in  the  vapor  pressure  drop  and  vapor  velocity  caused  by  the 
reduction  in  flow  area  would  not  have  a  significant  effect  on  heat  pipe 
capability. 
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CHAPltJ  5 


TRANSIENT  BEHAVIOR  OF  HP/TES  SYSTEM 
UNDER  REVERSED- PULSE  HEAT  LOADS  APPLIED  AT  THE  CONDENSER 


The  transient  behavior  of  a  high- temperature  axially  grooved  heat 
pipe  (HP)  which  incorporates  thermal  energy  storage  (TES) ,  under 
reversed- pulse  heat  loads  applied  at  the  condenser  is  presented  in  this 
chapter.  Liquid  sodium,  which  is  used  to  remove  the  heat  released  by  a 
power  generator,  circulates  through  a  HP/TES  cooling  device  from  which 
the  heat  is  rejected  into  space.  The  transient  response  of  three 
different  HP/TES  configurations  under  reversed- pulse  heat  loads  are 
compared:  (1)  a  heat  pipe  with  a  large  empty  cylinder  installed  in  the 
vapor  core,  (2)  a  heat  pipe  with  a  single  large  PCM  cylinder,  and  (3)  a 
heat  pipe  with  six  small  PCM  cylinders.  The  results  for  a  heat  pipe 
with  and  without  an  adiabatic  section  will  be  presented,  respectively. 

5.1  Description  of  The  Problem 

Future  space  missions  will  involve  thermal  transport  devices  with 
the  ability  to  handle  reversed- pulse  heat  loads.  Figure  5.1  shows  a 
schematic  diagram  of  the  cooling  system  for  a  power  generator.  A 
certain  amount  of  heat  is  continuously  being  released  by  the  power 
generator  and  removed  by  the  liquid  sodium  loop.  The  sodium  loop 
circulates  through  the  HP/TES  cooling  device,  where  the  heat  is  rejected 
into  space.  Under  normal  conditions,  the  system  is  operating  at  steady- 


93 


Figure  5.1  Schematic  diagram  of  the  cooling  system  for  a  power 
generator 
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state.  Suddenly,  an  ident  laser  strikes  the  condenser  section  of  the 
HP/TES  cooling  device.  Under  such  severe  operating  conditions,  the  heat 
released  by  the  power  generator  can  no  longer  be  removed  by  the  HP/TES 
cooling  device,  and  the  reversed- pulse  heat  load  caused  by  the  incident 
laser  will  be  reversely  transferred  back  into  the  liquid  sodium  loop. 
Incorporation  of  thermal  energy  storage  (TES)  into  heat  pipe  rejection 
systems  can  be  a  promising  method  to  mitigate  damage  from  reversed- pulse 
heat  loads.  The  transient  response  of  different  heat  pipes  (with  or 
without  phase- change  material)  unuer  reversed- pulse  heat  loads  will  be 
studied  in  this  chapter. 


5.2  Analytical  Model 


The  numerical  model  used  in  this  chanter  to  predict  the  transient 
response  of  the  HP/TES  cooling  system  has  already  been  described  in 
Chapters  3  and  4.  In  the  numerical  solutio**,  the  heat  pipe  evaporator 
wall  surface  temperature  was  assumed  equal  to  the  sodium  loop 
temperature  because  the  surface  heat  transfer  coefficient  is  very  high. 
The  liquid  sodium  loop  temperature  can  be  predicted  by  the  following 
equation: 

'^loop  -iF  =  «g  -  fle  (6^) 

Since  this  study  focused  only  on  a  heat  pipe  unit,  it  should  be 
noticed  that  and  in  Eq.  (67)  are  the  total  sodi  im  loop  heat 
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capacitance  and  the  heat  rate  released  from  the  power  generator, 
respectively,  divided  by  the  total  number  of  heat  pipes  in  the  system. 

The  heat  rate  released  from  the  power  generator  Q  is  always  positive 

and  remains  constant.  is  defined  as  the  heat  rate  transferred  from 
the  sodium  loop  to  the  heat  pipe,  and  is  evaluated  from  the  temperature 
gradient  in  the  heat  pipe  wall.  It  can  become  negative  when  a  reversed 
heat  load  is  applied.  If  liquid  sodium  with  velocity  1  m/s  is 
circulating  through  the  heat  pipe,  a  very  high  surface  heat  transfer 
coefficient  (about  5*10^  V/m^K)  can  be  obtained.  Since  the  heat 
transfer  coefficient  is  so  high,  we  will  assume  that  the  heat  pipe 
evaporator  wall  surface  temperature  is  equal  to  the  sodium  loop 
temperature. 

Ve  also  applied  a  simple  lumped- heat- capacity  model  to  predict  the 
transient  behavior  of  the  heat  pipe  without  a  PCM.  The  average  heat 
pipe  temperature  was  predicted  by  the  following  equation: 

%  tF  =  1e-  "c  (8«) 

where 

^ef'^hp  -  Irev* 

Eq.(68)  is  coupled  with  Eq.(67)  to  calculate  the  sodium  loop  and 
heat  pipe  temperatures.  For  the  lumped- heat- capacity  model,  an  average 
surface  heat  transfer  coefficient  between  the  liquid  sodium  loop  and  the 
heat  pipe  evaporator  of  h  =  5*10  w/ra  K  was  assumed.  The  heat 
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capacitance  of  the  heat  pipe  C 


hp 


is  about  540  J/K. 


5.3  Results  for  Heat  Pipes  With  Adiabatic  Section 

Ve  applied  the  numerical  model  to  a  grooved  heat  pipe  for  which  the 
specifications  were  given  in  Chapter  4.  The  total  length  of  the  heat 
pipe  was  1.0  m  with  an  adiabatic  section  having  length  of  0.4  m.  Under 
normal  condition,  heat  was  transferred  from  the  sodium  loop  to  the  heat 
pipe  evaporator  section  by  forced  convection  and  was  removed  at  the 
condenser  by  radiation. 

Figure  5.2  shows  the  transient  response  of  three  HP/TES 
configurations  with  =  1000  J/K  when  a  reversed  heat  load  is 

suddenly  applied  at  the  condenser.  Prior  to  t=  10  s,  the  three  heat 
pipes  all  are  operating  at  steady- state  conditions  with  the  temperature 
of  the  sodium  loop  maintained  at  950  K.  Under  this  steady- state 

condition,  the  total  heat  rate  transferred  from  the  sodium  loop,  Q^,  is 
about  0.78  kV  (for  an  average  surface  heat  flux  of  about  4.3  V/cm  )  and 
is  equal  to  the  total  heat  rate  removed  at  the  condenser  by  radiation 

heat  transfer,  Q  .  The  average  heat  pipe  temperature  is  about  945  K.  A 
heat  transfer  rate  of  0.78  kV  is  continuously  released  from  the  power 
generator  to  the  sodium  loop  throughout  the  entire  operating  period. 
After  t=  10  s,  a  reversed  heat  load  of  qj.gy=  10  V/cm  is  suddenly 
applied  to  the  condenser.  As  one  can  see,  the  temperature  of  the  heat 
pipe  without  PCM  increases  rapidly.  On  the  other  hand,  the  temperatures 
of  the  other  two  heat  pipes  with  PCM  also  increase  rapidly  immediately 


97 


TES 


98 


Figure  5.2  Transient  response  of  heat  pipe  with  Cjppp=  1000  J/K  under  a  reversed 
heat  load 


after  the  reversed  heat  load  is  applied.  But  the  rapid  temperature 
increase  is  arrested  after  the  PCM  reaches  its  melting  temperature  and 
starts  to  melt.  It  can  also  be  seen  that  the  temperature  increase 
during  the  melting  process  with  six  small  PCM  cylinders  is  slower  than 
with  one  large  PCM  cylinder.  This  slower  temperature  rise  occurs 
because  the  total  surface  area  of  the  six  small  PCM  cylinders  available 
for  heat  transfer  is  larger  and  the  heat  conduction  path  is  shorter 
compared  to  the  case  with  one  large  PCM  cylinder.  However,  the  six 
small  PCM  cylinders  will  be  completely  melted  earlier,  at  about  t=  80  s, 
compared  to  t=  110  s  for  the  one  large  PCM  cylinder.  After  the  PCM  is 
completely  melted,  both  PCM  equipped  heat  pipes  undergo  rapid 
temperature  increases  until  they  reach  a  new  steady- state  condition. 

Due  to  the  complexity  of  the  present  heat  pipe  problem,  we  are  not 
able  to  predict  the  error  of  the  numerical  solutions  accurately  because 
no  analytical  solution  is  available.  However,  in  order  to  retain  high 
accuracy,  the  time  step  and  grid  spacings  used  in  this  research  were 
chosen  to  satisfy  the  stability  and  accuracy  conditions  in  Ref.  [18]. 
Under  these  conditions,  the  error  of  the  present  numerical  solutions 
should  be  only  a  few  percent.  To  further  validate  the  numerical 
solutions,  we  applied  the  numerical  model  to  the  same  problem  depicted 
in  Fig  5.2  by  using  a  smaller  time  step  At=  0.025  s  and  reducing  the 
grid  spacings  in  the  r  and  z  directions  by  half  for  both  the  heat  pipe 
and  PCM.  As  we  can  see  from  Fig  5.2,  the  solutions  remain  almost  the 
same  with  smaller  values  of  the  time  step  and  grid  spacings. 

Figure  5.3  shows  the  variations  of  heat  input,  Q^,  and  heat  output, 
Qp,  for  three  different  HP/TES  configurations  with  CjQQp=  1000 
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J/K  under  a  reversed  heat  load  applied  at  the  condenser.  The  total  heat 

rate  transferred  from  the  sodium  loop  to  the  heai.  pipe  evaporator,  Q  , 
can  be  determined  by  knowing  the  temperature  gradient  inside  the  wall 
along  the  evaporator.  The  total  heat  rate  removed  from  the  condenser, 

Q  ,  is  the  summation  of  the  heat  removed  by  radiation  and  the  reversed 
heat  load.  Prior  to  t=  10  s,  all  three  heat  pipes  are  operating  at 

steady- state  with  a  heat  input  of  Q  =  0.78  kV  equal  to  the  heat  output 

*  2 
of  Q  =  0.78  kV.  After  a  reversed  heat  load  q_  =  10  V/cm  is  applied  at 

t=  10  s,  the  heat  outputs  at  the  condensers  of  the  three  heat  pipes  all 

become  negative,  indicating  that  there  are  external  heat  loads  being 

added  at  the  condensers.  However,  these  heat  outputs  all  begin  to 

increase  due  to  greater  heat  removal  by  radiation  at  the  higher 

condenser  wall  temperatures.  The  variation  of  heat  output  is  similar  to 

that  of  the  heat  pipe  temperature  because  the  heat  output  depends 

strongly  on  the  condenser  wall  surface  temperature. 

The  heat  input  variation  is  a  strong  function  of  the  sodium  loop 
heat  capacitance  After  the  reversed  heat  load  is  applied,  the 

heat  input  of  the  heat  pipe  without  PCI  decreases  very  rapidly  in  the 
first  10  seconds  and  is  reversed  to  negative.  After  t=  20  s,  the  heat 
input  begins  to  increase  because  the  heat  output  increases,  and  the 
reversed  heat  flow  effect  becomes  less  and  less.  The  heat  inputs  of  the 
other  two  heat  pipes  fitted  with  PCI  also  decrease  rapidly  immediately 
after  the  reversed  heat  loads  are  applied.  After  the  PCM  starts  to  melt 
at  about  t=  15  s,  the  heat  inputs  increase  very  rapidly.  This  is 
because  most  of  the  reversed  heat  load  is  absorbed  by  the  PCM  so  that 
the  heat  pipe  temperature  increase  becomes  very  slow.  However,  the  heat 
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inputs  decrease  slightly  during  later  stages  of  the  melting  process 
because  the  PCM  capability  to  absorb  the  reversed  heat  load  is 
declining.  The  heat  inputs  for  the  two  heat  pipes  equipped  with  PCM 
drop  again  and  reverse  to  negative  after  the  PCM  is  completely  melted 
because  the  heat  pipe  temperature  starts  to  increase  rapidly.  It  can 
also  be  seen  that  the  heat  input  of  the  heat  pipe  with  six  small  PCM 
cylinders  is  higher  than  that  of  the  one  with  a  single  large  PCM 
cylinder  during  the  melting  process.  The  greater  heat  input  occurs 
because  the  temperature  of  the  heat  pipe  with  six  small  PCM  cylinders  is 
lower.  The  three  different  HP/TES  configurations  all  tend  to  reach  a 
new  steady- state  temperature  with  same  original  heat  input  and  heat 
output  equal  to  0.78  kV. 

Figure  5.4  shows  the  transient  response  of  three  different  HP/TES 
configurations  with  10000  J/K  under  a  reversed  heat  load  suddenly 

applied  at  the  condenser.  Before  t=  10  s,  all  three  heat  pipes  are 
operating  at  steady- state  conditions  as  we  mentioned  in  the  earlier 
case.  After  t=  10  s,  a  reversed  heat  load  of  q_  =  10  V/cm  is  applied 
at  the  condenser.  As  can  be  seen  from  the  figure,  the  temperature 
increase  of  all  three  HP/TES  configurations  is  very  slow.  Because  of 
its  high  heat  capacitance,  the  sodium  loop  acts  like  a  huge  heat  sink 
which  can  absorb  most  of  the  reversed  heat  loads  and  arrest  the  heat 
pipe  temperature  increase.  It  is  clear  that  with  such  a  high  sodium 
loop  heat  capacitance,  installation  of  PCM  to  mitigate  the  reversed  heat 
loads  is  unnecessary. 

Figure  5.5  shows  the  variations  of  heat  input  and  heat  output  of 
three  different  HP/TES  configurations  with  10000  J/K  under  a 

reversed  heat  load  applied  at  the  condenser.  Compared  with  the  results 
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ariations  of  heat  input  and  heat  output  of  heat  pipes  with  C,  = 
10000  J/K  under  a  reversed  heat  load  ^ 


shown  in  Fig  5.3  for  the  case  with  Cj^QQp=  1000  J/K,  the  heat  inputs  of 
all  three  heat  pipes  decrease  very  rapidly  and  all  are  reversed  to 
negative  after  the  reversed  heat  load  is  applied.  With  such  a  high 
sodium  loop  heat  capacitance »  the  sodium  loop  itself  behaves  like  a 
massive  heat  sink  and  can  absorb  the  reversed  heat  loads  easily.  The 
heat  inputs  of  the  two  heat  pipes  with  PCM  are  reversed  less  than  the 
one  without  PCM  after  the  condenser  heat  loads  are  applied.  However, 
they  immediately  drop  again  after  the  PCM  is  completely  melted. 

Figure  5.6  presents  the  transient  response  of  the  heat  pipe  without 
PCM  under  reversed  heat  loads  as  predicted  by  the  lumped- heat- capacity 
model.  It  can  be  seen  that  the  results  from  the  lumped  model  have  very 
good  agreement  with  those  from  the  finite- difference  method.  The  lumped 
model  can  predict  the  average  heat  pipe  temperature  and  the  heat  flow 
input/output  at  the  evaporator  and  condenser  very  well  for  the  heat  pipe 
without  PCM.  The  heat  pipe  temperature  predicted  by  the  lumped  model  is 
about  10  K  lower  than  that  obtained  from  the  finite- difference  method 
throughout  the  time  period  of  interest.  This  discrepancy  arises  because 
the  heat  removed  from  the  condenser  by  radiation  is  overestimated  by 
using  the  average  heat  pipe  temperature  as  the  condenser  wall  surface 
temperature . 

Figure  5.7  shows  the  axial  variation  of  vapor  mass  flow  rate  for 
two  different  HP/TES  configurations  with  1000  J/I.  At  t=  10  s, 

both  heat  pipes  are  operating  at  steady- state  conditions  with  forward 
heat  loads  applied  at  the  evaporators.  All  the  vapor  mass  flow  rates 
are  positive  along  the  two  units.  The  vapor  mass  flow  rates  increase  in 
the  evaporator  section,  remain  almost  constant  in  the  adiabatic  section, 
and  then  decrease  in  the  condenser  section.  As  can  be  seen  from  Fig 
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5.5,  the  vapor  flow  of  the  heat  pipe  without  a  PCM  is  totally  reversed 
at  t=  60  s  since  both  heat  input  and  heat  output  are  negative.  For  the 
heat  pipe  with  six  small  PCM  cylinders,  the  vapor  flow  becomes  two 
separate  flows  with  opposite  directions  because  the  heat  input  is 
positive  while  the  heat  output  will  be  negative.  Evaporation  occurs  at 
both  evaporator  and  condenser  sections,  while  the  vapor  condenses  in  the 
adiabatic  section  and  on  the  outside  surfaces  of  the  PCM  containers.  As 
shown  in  Fig  5.7,  the  vapor  mass  flow  rate  is  positive  only  in  the 
evaporator  section  and  part  of  the  adiabatic  section.  It  is  negative 
over  the  remainder  of  the  heat  pipe.  One  should  also  note  that  the 
vapor  mass  flow  rate  at  the  adiabatic  section  of  the  heat  pipe  with  six 
small  PCM  cylinders  changes  more  rapidly  than  it  does  when  no  PCM  is 
present.  This  difference  occurs  because  a  considerable  amount  of  vapor 
condenses  on  the  outside  surfaces  of  the  PCM  containers  during  the 
melting  process. 

Figures  5.8a,b  show  the  axial  variation  of  vapor  pressure  and 
temperature  for  two  different  HP/TES  configurations  with  1000 

J/K.  At  t=  10  s,  both  heat  pipes  are  operating  at  steady- state 
conditions.  The  variation  of  the  vapor  pressure  for  both  heat  pipes  is 
almost  identical  since  there  is  little  difference  between  the  heat  pipe 
temperatures  and  vapor  mass  flow  rates.  At  t=  60  s,  both  heat  pipes 
have  a  higher  vapor  pressure  at  the  condenser  end  because  the  vapor 
flows  are  reversed.  The  vapor  pressure  drop  along  each  heat  pipe  at  t= 
60  s  is  much  less  than  at  t=  10  s.  For  the  heat  pipe  without  a  PCM,  the 
smaller  pressure  drop  is  mainly  due  to  its  higher  temperature,  which 
strongly  influences  vapor  pressure  variation  as  discussed  in  Chapter  4. 
For  the  heat  pipe  with  six  small  PCM  cylinders,  the  lower  pressure  drop 
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Figure  5.88  Axial  variation  of  v?nor  pressure  for  the  same  operating  o 
depicted  in  Fig.  5.‘/ 
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variation  of  vapor  temperature  for  the  same  operating 
ions  depicted  in  Fig.  5.2 


is  a  result  of  higher  heat  pipe  operating  temperature,  smaller  vapor 
mass  flow  rate  and  the  shorter  vapor  flow  path  resulting  from  the  two 
separate  vapor  flows.  One  should  also  note  that  the  axial  variation  of 
vapor  temperature  is  very  similar  to  that  of  vapor  pressure. 

The  transient  responses  of  the  heat  pipes  under  a  reversed- pulse 
heat  load  applied  to  the  condenser  from  t=  20  s  to  t=  80  s  are  shown  in 
Figs  5.9  and  5.10.  As  is  apparent  from  Fig  5.9,  the  temperature  of  all 
three  heat  pipes  responds  very  rapidly  and  starts  to  decrease  as  soon  as 
the  reversed- pulse  heat  load  is  removed  at  t=  80  s.  The  temperature  of 
the  heat  pipe  without  PCM  decreases  very  rapidly  after  the 
reversed- pulse  heat  load  is  removed.  The  temperature  of  each  unit  with 
PCM  also  decreases  rapidly  immediately  after  this  time,  but  the  decrease 
becomes  very  slow  when  the  PCM  reaches  the  melting  point  and  starts  to 
solidify.  The  six  small  PCM  cylinders  will  completely  solidify  earlier 
than  a  single  large  PCM  cylinder.  After  the  PCM  has  completely 
solidified,  the  temperature  of  both  heat  pipes  with  PCM  resumes  its 
decrease,  and  the  heat  pipes  gradually  return  to  their  initial 
steady- state  operating  conditions.  Approximately  1,660  seconds  are 
required  for  the  heat  pipe  with  six  small  PCM  cylinders  to  return  to  the 
initial  steady- state  condition  after  the  reversed- pulse  heat  load  is 
removed,  while  the  configuration  with  a  single  large  cylinder  needs 
2,260  seconds.  Fig  5.10  shows  the  percentage  of  PCM  melted  versus  time 
for  the  same  case  depicted  in  Fig  5.9.  One  can  readily  see  that  the  PCM 
does  not  respond  as  rapidly  as  the  heat  pipe  temperature  does.  In  fact, 
after  the  reversed- pulse  heat  load  is  removed  at  t=  80  s,  147.  more  of 
the  six  small  PCM  cylinders  and  257.  more  of  the  one  large  PCM  cylinder 
will  still  be  melted  before  the  phase  change  ceases.  It  is  also  clear 
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Figure  5,&  Transient  response  of  heat  pipes  under  a  reversed- pulse  heat  load 
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Figure  5. 1C  Portion  of  PCM  melted  versus  time  for  heat  pipes  under 
reversed- pulse  heat  load 


that  the  small  PCM  cylinders  solidify  faster  than  the  larger  one  does. 

Figures  5.11  and  5.12  illustrate  the  results  for  the  transient 
response  of  the  heat  pipes  with  periodic,  reversed- pulse  heat  loads. 

The  time  period  is  2000  s  and  each  of  the  reversed- pulse  heat  loads 
lasts  60  s.  The  temperature  response  in  each  time  period  is  similar  to 
the  results  shown  in  Fig  5.9.  The  temperature  of  the  heat  pipe  without 
a  PCM  simply  oscillates  up  and  down  periodically.  The  temperature  of 
the  heat  pipe  with  six  small  PCM  cylinders  remains  almost  constant 
throughout  the  whole  period  due  to  the  very  efficient  melting  and 
solidification  of  this  PCM  configuration.  As  shown  in  Fig  5.12,  the 
molten  fraction  of  the  PCM  for  the  heat  pipe  with  six  small  PCM 
cylinders  does  not  differ  from  one  cycle  to  the  next.  This  is  because 
the  six  small  PCM  cylinders  solidify  completely  at  the  end  of  each 
cycle.  The  percentage  of  PCM  melted  for  the  heat  pipe  with  one  large 
PCM  cylinder  continues  to  increase  as  the  pulse  cycles  continue.  This 
occurs  because  the  one  large  PCM  cylinder  does  not  have  sufficient  time 
to  solidify  completely  at  the  end  of  each  time  period  for  this  cyclic 
rate. 

Figures  5.13  and  5.14  present  the  results  for  the  heat  pipes  under 
periodic,  reversed- pulse  heat  loads  with  a  shorter  time  period  of  600  s. 
Each  of  the  reversed- pulse  heat  loads  also  lasts  60  s.  It  is  clear  from 
these  figures  that  the  temperature  of  both  heat  pipes  with  PCM  is  under 
control  during  the  first  cycle.  However,  the  temperatures  become  very 
high  during  the  next  two  cycles.  This  failure  of  the  PCM  temperature 
regulation  mechanism  occurs  because  the  time  period  is  too  short  for  the 
PCM  to  solidify  completely  at  the  end  of  each  cycle.  At  the  end  of  the 
first  cycle,  only  about  507.  of  the  PCM  is  solidified  for 
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Figure  5.12  Portion  of  PCM  aelted  versus  tine  for  heat  pipes  under  a  periodic 
reversed- pulse  heat  load  with  tine  period  of  2000  s 
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Figure  5.14  Portion  of  PCM  aelted  versus  time  for  heat  pipes  under  a  periodic 
reversed- pulse  heat  load  with  time  period  of  600  s 


either  heat  pipe.  After  the  second  reversed- pulse  heat  load  is  applied, 
the  PCM  soon  becomes  completely  melted,  and  the  heat  pipe  temperature 
increases  rapidly.  The  PCI  will  never  again  have  time  to  solidify  and 
thus  becomes  useless  for  thermal  regulation  purposes  after  the  second 
cycle.  Obviously,  cycle  time  is  a  vital  consideration  in  the  design  of 
a  HP/TES  cooling  system  to  handle  periodic,  reversed- pulse  heat  loads. 

5.4  Results  for  Heat  Pines  Without  Adiabatic  Section 

In  this  section,  the  numerical  model  will  be  applied  to  a  grooved 
heat  pipe  without  an  adiabatic  section  under  partially  reversed  heat 
loads  applied  at  the  condenser.  The  specifications  of  the  heat  pipe  are 
the  same  as  those  given  in  Chapter  4,  except  there  is  no  adiabatic 
section.  The  total  length  of  the  heat  pipe  is  1.0  m  with  the  evaporator 
and  condenser  sections  having  lengths  of  0.3  m  and  0.7  m,  respectively. 
The  numerical  model  predicted  the  transient  response  of  the  HP/TES 
cooling  system  under  a  variety  of  partially  reversed  heat  loads  as  shown 
in  Fig  5.15. 

Figure  5.16  shows  the  transient  response  of  three  HP/TES 
configurations  with  =  1000  J/K  when  a  reversed  heat  load  is 

suddenly  applied  to  757i  of  the  condenser  surface.  Prior  to  t=  10  s,  the 
three  heat  pipes  all  are  operating  at  steady-  state  conditions  with  the 
temperature  of  the  sodium  loop  maintained  at  950  K.  Under  this 
steady- state  condition,  the  total  heat  rate  transferred  from  the  sodium 

loop,  Qg,  is  about  1.78  kV  (for  an  average  surface  heat  flux  of  about 
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Figure  5.15  A  heat  pipe  without  an  adiabatic  section  operates  under  a 
partially  reversed  heat  load 
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Figurft  5. It  Transient  response  of  heat  pipes  under  a  reversed 
covers  757.  of  condenser 


2 

9.8  V/cm  )  and  is  equal  to  the  total  heat  rate  removed  at  the  condenser 

by  radiation  heat  transfer  Q^.  The  average  heat  pipe  temperature  is 
about  938  K.  A  heat  transfer  rate  of  1.78  kV  is  continuously  released 
from  the  power  generator  to  the  sodium  loop  throughout  the  entire 
operating  period.  After  t=  10  s,  a  reversed  heat  load  of  qj.gy=  20 
v/cm  is  suddenly  applied  to  757,  of  the  condenser  surface.  Compared 
with  the  results  shown  in  Fig  5.2,  the  PCM  is  completely  melted  earlier 
due  to  the  higher  heat  loads  applied  at  the  evaporator  and  condenser. 

The  six  small  PCM  cylinders  will  be  completely  melted  at  about  t=  40  s 
and  the  single  large  PCM  cylinder  at  t=  60  s.  As  one  can  see  from  the 
figure,  the  temperature  of  both  heat  pipes  fitted  with  PCM  also 
increases  very  swiftly  during  the  PCM  melting  process.  This  rapid 
temperature  rise  happens  because  that  the  radial  temperature  gradient 
inside  the  PCM  is  very  large  under  such  high  heat  loads.  The  surface 
temperature  of  a  single  large  PCM  cylinder  is  about  90  K  higher  than  the 
PCM  melting  point  at  t=  50  s.  In  designing  an  HP/TES  cooling  system  to 
handle  such  high  total  reversed  heat  loads,  one  must  install  more  PCM 
cylinders  in  a  heat  pipe  with  a  larger  vapor  flow  area.  Thus  the  heat 
load  to  each  PCM  cylinder  can  be  reduced,  the  radial  temperature 
gradient  inside  the  PCM  can  be  decreased,  and  the  PCM  melting  process 

will  last  longer.  The  variations  of  the  heat  input,  Q^,  and  heat 

output,  Q  ,  as  shown  in  Fig  5.17  are  similar  to  the  results  already 
given  in  Fig  5.3. 

Figures  5.18  and  5.19  present  the  transient  response  of  three 
HP/TES  configurations  with  =  1000  J/K  when  a  reversed  heat  load  is 

suddenly  applied  to  507,  of  the  condenser  surface.  The  six  small  PCM 
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Figure  5.17  Variations  of  heat  input  and  heat  output  of  heat  pipes  under 
reversed  heat  load  which  covers  757.  of  condenser 
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Figure  5.18  Transient  response  of  heat  pipes  under  a  revcr 
covers  507i  of  condenser 


Figure  5.16  Variations  of  heat  input  and  heat  output  of  heat  pipes  under- 
reversed  heat  load  which  covers  507,  of  condenser 


cylinders  and  the  one  larger  PCM  cylinder  will  be  completely  melted  at 
t=  50  s  and  80  s,  respectively.  Although  both  the  fraction  of  the 
condenser  surface  directly  exposed  to  the  reversed  heat  loads  and  the 
total  load  itself  are  reduced  from  the  previous  case,  the  temperature  of 
both  heat  pipes  still  increases  rapidly  during  the  PCM  melting  process. 
The  surface  temperature  of  the  single  large  PCM  cylinder  is  about  80  K 
higher  than  the  PCM  melting  point  at  t=  70  s. 

Figures  5.20  and  5.21  show  the  transient  response  of  three  HP/TES 
configurations  with  =  1000  J/K  when  a  reversed  heat  load  is 

suddenly  applied  to  257,  of  the  condenser  surface.  The  six  small  PCM 
cylinders  and  the  single  large  PCM  cylinder  will  now  be  completely 
melted  at  t=  80  s  and  120  s,  respectively.  As  one  can  see,  the 
temperature  of  the  heat  pipe  with  six  small  PCM  cylinders  only  increases 
about  10  K  during  the  PCM  melting  process.  The  surface  temperature  of 
the  single  large  PCM  cylinder  is  about  40  K  higher  than  the  PCM  melting 
point  at  t=  110  s;  however,  the  surface  temperature  excess  is  only  about 
12  K  for  the  small  PCM  cylinders  at  t=  70  s. 

Figure  5.22  shows  the  axial  variation  of  vapor  mass  flow  rate  for 
the  same  case  as  shown  in  Fig  21.  At  t=  10  s,  both  heat  pipes  are 
operating  at  steady- state  conditions  with  forward  heat  loads  applied  at 
the  evaporators.  All  the  vapor  mass  flow  rates  are  positive  along  the 
two  units.  The  vapor  mass  flow  rate  increases  at  the  evaporator  section 
and  then  decreases  at  the  condenser  section.  At  t=  60  s,  each  vapor 
flow  in  both  heat  pipes  becomes  two  separate  flows  moving  in  opposite 
directions.  Evaporation  occurs  at  both  the  evaporator  and  the  portion 
of  the  condenser  section  where  the  reversed  heat  load  is  applied,  while 
the  vapor  condenses  in  the  remaining  portion  of  the  condenser  section 
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Transient  response  of  beat  pipes  under  a  reversed  heat  load  which 
covers  2575  of  condenser 


128 


Figure  5.21  Variations  of  heat  input  and  heat  output  of  heat  pipes  under 
reversed  heat  load  which  covers  257.  of  condenser 


covers  257.  of  condenser 


and  on  the  outside  surfaces  of  the  PCM  containers. 

The  transient  responses  of  the  heat  pipes  under  a  reversed- pulse 
heat  load  applied  to  257i  of  the  condenser  surface  from  t=  20  s  to  t=  80 
s  are  shown  in  Figs  5.23  and  5.24.  The  trends  depicted  by  the  transient 
responses  are  similar  to  the  results  shown  in  Figs  5.9  and  5.10. 

However,  compared  to  the  results  in  Figs  5.9  and  5.10,  the  PCM  cylinders 
are  completely  solidified  much  more  quickly  after  the  reversed- pulse 
heat  loads  are  removed.  Only  about  1,000  seconds  are  required  for  the 
heat  pipe  with  six  PCM  cylinders  to  return  to  its  initial  steady- state 
condition  after  the  reversed- pulse  heat  load  is  removed,  and  the 
configuration  with  one  large  cylinder  can  accomplish  this  in  about  1,270 
seconds.  This  shorter  adjustment  time  can  be  achieved  because  the  heat 
will  be  removed  more  quickly  from  a  heat  pipe  with  a  larger  condenser 
area. 
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which  covers  257*  of  condenser 
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Figure  5.24  Portion  of  PCM  aelted  versus  time  for  heat  pipes  under  a 
reversed- pulse  heat  load  which  covers  257,  of  condenser 


CHAPTEE  6 


CONCLUSIONS 

In  this  research,  the  transient  behavior  of  a  high- temperature 
grooved  heat  pipe  with  thermal  energy  storage  (TES)  under  pulse  heat 
loads  was  modeled  using  a  three-dimensional  ADI  finite- difference 
method.  A  phase- change  material  (PCM)  encapsulated  in  cylindrical 
containers  was  used  as  the  TES.  Two  different  types  of  pulse  heat  loads 
applied  to  the  heat  pipe  were  studied.  One  was  high- pulse  heat  loads 
applied  at  the  heat  pipe  evaporator.  The  other  one  was  reversed- pulse 
heat  loads  applied  at  the  condenser.  It  was  found  that  the  PCM  is  very 
effective  in  mitigating  the  adverse  effects  of  both  types  of  pulse  heat 
loads.  The  six  small  PCM  cylinders  are  more  efficient  than  the  single 
large  PCM  cylinder  in  reducing  the  rapid  increase  in  heat  pipe 
temperature  under  pulse  heat  loads,  and  they  can  also  handle  periodic, 
pulse  heat  loads  better  since  they  solidify  faster. 

Since  the  heat  pipe  capillary  limit  is  dependent  on  the  overall 
liquid  and  vapor  pressure  drops  and  the  vapor  pressure  drop  dominates 
the  liquid  pressure  drop  for  this  type  of  heat  pipe,  the  effect  of  the 
liquid  pressure  drop  on  the  heat  pipe  capillary  limit  can  be  neglected. 
The  vapor  pressure  and  temperature  drops  of  the  heat  pipe  were  found  to 
be  strongly  dependent  on  the  operating  vapor  temperature.  The  vapor 
flow  can  be  reversed  or  become  two  opposing  flows  under  reversed  heat 
loads.  The  numerical  results  also  indicated  that  the  heat  inputs  and 
outputs  of  the  heat  pipes  in  the  cooling  system  under  reversed- pulse 
heat  loads  are  strongly  dependent  on  the  sodium  loop  heat  capacitance. 

The  main  disadvantage  of  installing  phase- change  material  in  the 
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vapor  core  of  a  heat  pipe  is  the  accompanying  reduction  in  the  vapor 
flow  area.  This  reduction  in  vapor  flow  area  could  cause  the  vapor 
pressure  drop  and  vapor  velocity  to  increase  significantly,  and  thus 
decrease  the  heat  pipe  capability.  Fortunately,  the  PCM  itself  can 
absorb  a  large  portion  of  the  heat  loads  during  the  melting  process 
after  pulse  heat  loads  are  applied.  Some  vapor  will  condense  on  the 
surface  of  PCM  cylinders  and  reduce  the  vapor  velocity  and  pressure 
drop.  From  the  numerical  results,  one  can  find  that  not  only  can  the 
PCM  compensate  for  the  decrease  in  heat  pipe  transport  capability  due  to 
the  reduction  in  vapor  flow  area,  it  can  also  actually  handle  the  pulse 
heat  loads  very  effectively. 

In  the  design  of  an  HP/TES  system,  one  should  choose  a  PCM  with  a 
melting  point  slightly  higher  than  the  normal  operating  temperature. 

Then  if  a  pulse  heat  load  higher  than  the  heat  pipe  transport  limitation 
is  applied,  the  PCM  can  respond  fast  enough  to  begin  melting  and  absorb 
some  of  the  heat  before  the  heat  pipe  reaches  its  operating  limit  and 
bums  out.  To  reduce  the  chance  of  complete  melting  during  the  pulse 
period,  the  latent  heat  of  fusion  of  the  PCM  should  be  as  large  as 
possible.  Also,  the  cycle  time  is  a  vital  consideration  in  the  design 
of  an  HP/TES  cooling  system  to  handle  periodic,  pulse  heat  loads. 

The  concept  of  incorporating  phase- change  material  inside  a 
low- temperature  heat  pipe  (such  as  a  water  heat  pipe)  is  also  sound  if 
the  goal  is  to  limit  the  temperature  extremes  encountered  when  the  heat 
load  is  time  dependent.  For  most  low- temperature  heat  pipes,  the  vapor 
pressure  drop  is  small,  and  vapor  flow  usually  does  not  play  an 
important  role  in  determining  the  heat  pipe  transport  capability.  Thus 
the  increase  in  the  vapor  pressure  drop  and  vapor  velocity  caused  by  the 
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reduction  in  flow  area  would  not  have  a  significant  effect  on  heat  pipe 
capability. 


An  improved  three-dimensional  ADI  finite- difference  method  used  to 
model  the  heat  conduction  through  the  heat  pipe  wall  and  wick  was 
developed  in  this  research.  An  important  characteristic  of  this  new  ADI 
method  is  that  the  resulting  finite- difference  equations  are  consistent 
with  physical  considerations.  Compared  to  the  conventional  ADI  method, 
this  modification  allows  the  time  step  to  be  increased  by  about  2  orders 
of  magnitude  without  compromising  significantly  on  the  accuracy  of  the 
numerical  solution.  Compared  with  the  well-known  Brian  and  Douglas  ADI 
methods,  this  new  ADI  method  yields  higher  accuracy  and  requires  less 
computer  storage. 

The  equivalent  heat  capacity  method  proposed  by  Hsiao  for  the 
Stefan  problem  was  tested  in  this  research,  but  a  large  energy  balance 
error  was  found.  Similar  conclusion  was  also  given  by  Pham.  Comparing 
with  exact  solutions,  Pham  pointed  out  that  the  Hsiao  method  yields 
results  with  up  to  227,  error.  The  low  accuracy  of  the  Hsiao’s  method 
could  be  due  to  its  ambiguous  theoretical  basis.  Pham  suggested  a 
simple  and  accurate  method  which  includes  the  good  features  of  the 
enthalpy  and  heat  capacity  methods.  One  of  the  good  features  in  the 
Pham  method  is  that  it  estimates  the  new  temperature  from  the  estimated 
enthalpy  change  to  qvoid  the  problem  of  jumping  the  latent  heat  peak. 
However,  the  Pham  method  has  a  singularity  problem. 

In  this  report,  we  adopted  the  best  features  of  the  Pham  method  for 
the  Stefan  problem  and  made  some  modifications  to  improve  on  its  weak 
points.  This  modified  Pham  method  was  used  in  conjunction  with  a 
two-dimensional  ADI  scheme.  Compared  with  analytical  solutions,  the 
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present  method  for  melting  and  solidification  was  found  to  have  very 
good  accuracy  without  the  singularity  problem  of  the  Pham  method. 

Liquid  dynamics  in  a  heat  pipe  become  very  significant  when  dryout 
and  rewetting  occur  in  the  heat  pipe  wick.  However,  dryout  depends  on 
the  heated  zone  and  the  instantaneous  local  saturation.  Thus,  if  dryout 
is  to  be  accurately  predicted,  the  temporal  dependence  of  the  saturation 
distribution  must  be  taken  account.  A  complete  capillary  liquid  flow 
model  is  needed  to  predict  the  dryout  and  rewetting  behaviors  of  a  heat 
pipe.  A  good  capillary  liquid  flow  model  should  include  the  effect  of 
vapor  pressure  changes  on  the  liquid  meniscus  contact  angle  and  also  be 
able  to  predict  the  saturation  distribution. 

It  was  found  that  the  lumped- heat- capacity  model  can  predict  the 
average  heat  pipe  temperature  and  the  heat  flow  input/output  at  the 
evaporator  and  condenser  very  well  for  the  heat  pipe  without  PCM.  Ve 
have  also  checked  the  total  energy  balance  and  found  the  error  is  less 
than  17i  at  each  time  step.  To  further  validate  the  numerical  solutions, 
we  applied  the  numerical  model  to  some  problems  by  using  different  time 
steps  and  grid  spacings.  The  results  reported  here  were  obtained  with  a 
sufficiently  fine  grid  spacing  and  time  step  so  that  numerical  results 
are  essentially  independent  of  the  time  step  and  grid  spacing. 
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APPENDH 

Listing  of  The  Computer  Program 
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C  3-D  COMPUTER  PROGRAM  FOR  TRANSIENT  BEHAVIOR  OF  AXIALLY  GROOVED  HEAT 
C  PIPE  UNDER  PULSE  HEAT  LOADS 
C 
C 

C  CKL  -THERMAL  CONDUCTIVITY  OF  LIQUID  (W/M-K) 

C  CKLHL  -THERMAL  CONDUCTIVITY  OF  LIQUID  LIH  (W/M-K) 

C  CKLHS  -THERMAL  CONDUCTIVITY  OF  SOLID  LIH  (W/M-K) 

C  CKSS  -THERMAL  CONDUCTIVITY  OF  SOLID  'W/M-K) 

C  CKP  -NODAL  THERMAL  CONDUCTIVITY  (W/M-K) 

C  CL  -LIQUID  HEAT  CAPACITY  (J/K-KG) 

C  CLHL  -HEAT  CAPACITY  OF  LIQUID  LIH  (J/K-KG) 

C  CLHS  -HEAT  CAPACITY  OF  SOLID  LIH  (J/K-KG) 

C  CP  -NODAL  HEAT  CAPACITY  (J/K-KG) 

C  CS  -SOLID  HEAT  CAPACITY  (J/K-KG) 

C  DENL  -DENSITY  OF  LIQUID  (KG/M3) 

C  DENLHL  -DENSITY  OF  LIQUID  LIH  (KG/M3) 

C  DENLHS  -DENSITY  OF  SOLID  LIH  (KG/M3) 

C  DENP  -NODAL  DENSITY  (KG/M3) 

C  DENS  -DENSITY  OF  SOLID  (KG/M3) 

C  DQVDA(K)  -LOCAL  INWARD  HEAT  FLUX  ON  THE  VAPOR  INTERFACE  (W/M2) 

C  DRS  -SPACE  INCREMENT  IN  RADIAL  DIRECTION  IN  SOLID  (M 

C  DRW  -SPACE  INCREMENT  IN  RADIAL  DIRECTION  IN  WICK  (M) 

C  DT  -TIME  STEP  (SEC) 

C  DTHETA  -ANGULAR  INCREMENT  (RADIANS) 

C  DTVl  -HALF  TEMPERATURE  INTERVAL  OF  VAPOR  FOR  SEARCHING  TV(1)  (K) 

C  D2  -SPACE  INCREMENT  IN  AXIAL  DIRECTION  (M) 

C  EMISS  -EMISSIVITY 

C  HFUSLH  -FUSION  HEAT  OF  LIH  (J/KG) 

C  HTCAMB  -AMBIENT  HEAT  TRANSFER  COEFFICIENT  FOR  THE  CONDENSER  (W/M2K) 

C  NG  -NUMBER  OF  GROOVES 

C  PV(K)  -VAPOR  THERMDYNAMIC  PRESSURE  (N/M2) 

C  QEVAP  -HEAT  FLUX  ON  THE  WALL  OF  THE  EVAPORATOR  (W/M2) 

C  QVIN  -NODAL  HEAT  TRANSFER  FROM  HEAT  PIPE  WALL  TO  VAPOR  FLOW  (W) 

C  QV(K)  -LOCAL  HEAT  TRANSFER  FROM  HEAT  PIPE  WALL  TO  VAPOR  FLOW  (W) 

C  R  -NODAL  RADIUS  (M) 

C  RI  -INSIDE  RADIUS  OF  HEAT  PIPE  WALL  (M) 

C  RMENIS  -MENISCUS  RADIUS  (M) 

C  RO  -OUTSIDE  RADIUS  OF  HEAT  PIPE  WALL(M) 

C  RW  -VAPOR  CORE  RADIUS  (M) 

C  SKsMA  -STEFAN-BOLTZMAN  CONSTANT  (W/M2K4) 

C  SUMQVI  -TOTAL  EVAPORATION  RATE  (W) 

C  SUMQVO  -TOTAL  CONDENSATION  RATE  (W) 

C  TAMB  -AMBIENT  TEMPERATURE  TO  THE  CONDENSER  (K) 

C  TMELLH  -MELT  TEMPERATURE  OF  LIH  (K) 

C  T(I,J,K)  -NODAL  TEMPERATURE  (K) 

C  TV(K)  -VAPOR  TEMPERATURE  (X) 

C  XL  -HEAT  PIPE  LENGTH  (M) 

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

REAL  *0  TVOLD(100),TV(100),TE(40,100),TTE(40,100),DTV(100) 

REAL  *8  T(8,5,100),TW(5,100),TT(8,5,100),TOLD(8,5,100) 

REAL  *8  TVO  (1(j0)  ,TVNEW{100)  ,DPV(100)  ,DPL(100) 
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o  o  u  u 


DIMENSION  P (100) ,0(100) ,CKNR(5, 100) 

DIMENSION  QWAL1,(5,100)  ,DQVDAN(100)  ,DQDATN(100) 
DIMENSION  THPAVG (20000) ,TVAVG (20000) ,CT1ME (20000) 
DIMENSION  QVO  (100)  ,QVNEW(100)  ,QVOIiD(100) 


COMMON  /HPIPE/  TV.DQDATEdOO)  ,  ZL,D2,NZ,NDT,QVTES  (100)  , 
SNTESTU, TIME, KBURN, NWRITE, XWRITE, NGROV, RW, NADIA, PV ( 100 ) 

COMMON  /VCORE/  DQVDA (100) , DQDABAR (100) , IPCM, PVTTL (100) 

$  ,TVREAL(100) ,PLM,DTES,DVAPOR,DWlCK,DHYDRO,XMMAX,KPLM 
$  ,NEVAP,XTI,DENV(100) 

COMMON  /PCM/  TVOLD,TE,TTE,DTV,DT,NADI,NRTES,DTHETA,RTES, 

$  TSTART,NTHETA,QPCM,QMEI,T,QPCMTTL 

COMMON  /LIQ/  QV(100),PL(100),SUMQVI,KPVMIN,PC»1AX,DPVPLMAX, 

$  CL, NZPl , WIDBAR, DEPTHG, DENL, VISLIQ, TENLIQ, LAMBDA 


CHOOSE  RUNNING  CASE:  IPCM-1  (WITHOUT  PO!) ;  IPCM-2  (WITH  PCM) 

METHOD-1  (EXPLICIT) ;  METHOD-2  (IMPLICIT) 


IPCM-1 

METHOD-1 


FRELAX-0. 02 

IF (METHOD. EQ.l)  FRELAX-0. 0 


C  CHOOSE  TIME  STEP  (SEC) 

DT-0 . 1 

NDTEND-  6000 

C  SPECIFY  THERMAL  ENERAGY  STORAGE 
NTESTU-  0 
RTES-0.004 
NRTES-  40 


CC  NTESTU-  6 

CC  RTES-0. 001633001 

CC  NRTES-  10 

NWINCH-200 


C  CH<X)SE  TIME  TO  WRITE 
NWRITE-  100 
XWRlTE-1 . 0*NWRITE 

C 
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CONVERGENCE  CRITERIA  FOR  THE  EVAPORATOR  END  VAPOR  TEMPERATURE 
EPSQVK-0 . 01 
C 
C 

C  SPECIFIED  TRUNCATION  ERROR  FOR  THE  VAPOR  FLOW 
EPSV-0.01 

C  SPECIFIED  BOUNDARY  CONDITIONS  OF  HEAT  TRANSFER  ON  THE  HEAT  PIPE  WALL 
QEVAP-  4.3*10000.0 
CC  TF-941.0 

CC  HTCF-37043.0 

EMISS-1.0 
SIGMA-5. 669E-8 
TAMB-  0 . 0 

KBURN-0 

C 

C 

C  HEAT  PIPE  GEOMETRY 

NGROV-18 
NTHETA-4 
NTHEPG-4 
NRS-4 
NRW-4 

NR-NRS+NRW 

NZ-  40 
NEVAP-12 
NADIA-  28 
ZL-  1.0 
NZPl-NZ+1 

PLM-NEVAP*DZ 

KPLM-NEVAP 


RO-0.009525 
RI-0. 008262 
RW-0. 007000 

NWIRE-NWINCH/0 . 0254 
RC-1.0/(2.0*NWIRE) 


WIDBAR-2 . 0*3 . 1416* (RI+RW) /2 . 0/2 . O/NGROV/2 . 0 

DEPTHG-RI-RW 

DWICK-2 . 0*RW 

DTES-2.0*RTES 


DRS-(RO-RI) /NRS 
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DRW- (RI-RW) /NRW 
DTHETA-2*3 . 14/NTHETA/NGROV 
DZ-ZL/NZ 

XTI-NEVAP*DZ 

C  GROOVE  GEOMETRY 
JING2-2 
J1NG3-3 

C  THERMAL  CONDUCTIVITY,  HEAT  CAPACITY  AND  DENSITY 
C  SOLID: STAINLESS  STEEL,  LIQUID: SODIUM 
CKSS-18.4 
CKL-75.7 
CS-530.0 
CL-1305.0 
DENS-7744.0 
DENL-882.0 
NTHEPl-NTHETA 


WRI TE ( 6 , * ) • NTESTU- * , NTESTU 
WRITE (6,*) 'RO-',RO, •PL-',ZL 

WRITE  (6,*)  'NZ-',NZ,  'NR-', NR,  'NRS-^NRS,  'NRW-', NRW 
WRITE(6,*) 'QEVAP-',QEVAP, • SIGMA-' , SIGMA, •TAMB-',TAMB 


C  HEAT  PIPE  STARTING  TEMPERATURE 
TSTART-  940.0 
DO  30  K-1,NZ 
DO  33  J-1,NTHEP1 
DO  36  I-1,NR 
T (I, J,K) -TSTART 
TOLD ( I , J, K) -TSTART 
TW(J,K) -TSTART 
36  CONTINUE 
33  CONTINUE 

TV (K) -TSTART 
TVOLD (K) -TSTART 
QV(K)-0.0 
DQVDA(K)-0.0 
QVTES(K)-0.0 
DQDATE (K) -0 . 0 
DTV(K)-0.0 
30  CONTINUE 

QETTL-0 . 0 
QCTTL-0 . 0 
QTES-0 . 0 
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SUMQVI-0.0 


TlME-0.0 

DO  60  NDT-  1,NDTEND 

C  PULSE  HEAT  FLUX  SPECIFICATION 
IF (NDT. LE. 100)  PULSE-1.0 
IF (NDT. GT. 100)  PULSE-10.0/4.3 


CC  IF (NDT. LE. 100)  PULSE-1.0 

CC  IF(NDT.GT.IOO.AND.NDT.LE.IOOO)  PULSE-10.0/4.3 

CC  IF(NDT.GT.IOOO)  PULSE-1.0 


CC  IF (NDT. LE. 200)  PULSE-1.0 

CC  IF(NDT.GT.200.AND.NDT.LE.400)  PULSE-10.0/4.3 

CC  IF(NDT.GT. 400. AND. NDT. LE. 2200)  PULSE-1.0 

CC  IF(NDT.GT. 2200. AND. NDT. LE. 2400)  PULSE-10.0/4.3 

CC  IF (NDT. GT. 24 00. AND. NDT. LE. 4200)  PULSE-1.0 

CC  IF (NDT. GT. 4200. AND. NDT. LE. 4400)  PULSE-10.0/4.3 

CC  IF (NDT. GT. 4400. AND. NDT. LE. 6200)  PULSE-1.0 

CC  IF(NDT.GT, 6200. AND. NDT. LE. 6400)  PULSE-10.0/4.3 

CC  IF(NDT.GT. 6400. AND. NDT. LE. 8200)  PULSE-1.0 

CC  IF(NDT.GT.8200.AND.NDT.LE.8400)  PULSE-10.0/4.3 

CC  IF(NDT.GT. 8400. AND. NDT. LE. 10200)  PULSE-1.0 

IF(NTESTU.EQ.O)  DHYDRO-0.014 
IF(NTESTU.EQ.l)  DHYDRO-0 . 0082 
IF(NTESTU.E(3.6)  DHYDRO-0 . 0075 

DVAPOR-(DWICK**2.0-NTESTU*DTES**2.0) **0.5 

IF (NDT. LE. 20)  EPSV-0.01 
IF (NDT. GT. 20)  EPSV-0.01 


F-0.01 

Fl-3.0-2.0*F 
C  DDTVl-1 . 0 

C  DTV1-(TV(1) -TVOLD(l) ) +DDTV1 

DTVl-3.0 

DO  10  K-1,NZ 
TVOLD(K)-TV(K) 

10  CONTINUE 

IF ( (NDT/XWRITE) .EQ. (NDT/NWRITE) )  THEN 


142 


WRITE(6,40)DT,F 

40  FORMAT (3X, 'THIS  OUTPUT  IS  FOR  DT-' , F6 . 3, 3X, ' F-' , F6 . 3// ) 
ENDIF 


TIME-NDT*DT 
C  LIQUID  PROPERTIES 

VISLlQ-0 . 0011-1 . 45E-6*TV (1) +5 . 32E-10*TV (1 ) **2 . 0 
TENLlQ-0 . 22-9 . 1E-5*TV (1) 

IiAMBDA-4636437 . 0-180 . 0*TV(1) 

PCMAX-2 . 0*TENLIQ/RC 

DO  20  K=1,NZ 
QVO (K) -QV (K) +QVTES (K) 

TV0(K)-TV(K) 

20  CONTINUE 

DO  25  K-1,NZ 

IF(K.LE.NEVAP)  DQDABAR (K) - (SUMQVI/NEVAP) / (2 . 0*3 . 1416*RW*DZ 
$  +NTESTU*2.0*3.1416*RTES*DZ) 

IF (K . GT , NEVAP . AND . K . LE . NADIA)  DQDABAR (K) -0 . 0 

IF  (K.GT, NADIA)  DQDABAR (K) —  ( SUMQVI /  (NZ-NADIA) )  /  (2 . 0*3 . 1416*RW*DZ 
$  +NTESTU*2.0*3.1416*RTES*DZ) 

25  CONTINUE 
C 

C  &&&&&&&  COMPUTATION  FOR  THE  FIRST  ONE-THIRD  TIME  STEP  &&&&&&&&& 
NADI-123 
QHP-0 . 0 

DO  180  J-1,NTHEP1 
NGORD- ( J-1 ) /NTHEPG+1 
JING-J- (NGORD-1) *NTHEPG 
DO  160  I«1,NR 
DO  100  K-1,NZ 

IF(NDT,EQ.l)  QV(K)-0.0 

IF(K.LE.KBURN)  THEN 
CKL-0. 0000001 
ELSE 

CKL-75,7 

ENDIF 


IF (K. LE. NEVAP)  QWALL (J, K) -PULSE*QEVAP 
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CC  IF(K.LE.NEVAP)  QWALL ( J, K) -HTCF* (TF-TW ( J, K) ) 

IF (K . GT . NEVAP . AND . K . LE . NADIA)  QWALL ( J, K) -0 . 0 

IF (K.GT. NADIA)  QWALL ( J, K) -EMISS*SIGMA* (TAMB**4 . 0-TW ( J, K) **4 . 0) 


IF(I.LE.NRS)  THEN 
R-RO- (1-0.5) *DRS 
DR-DRS 
CKP-CKSS 
CKN-CKSS 
CKS-CKSS 
CKE-CKSS 
CKW-CKSS 
CKF=CKSS 
CKB=CKSS 
DENP-DENS 
CP-CS 

I F ( I . EQ . NRS . AND . JING . GE . JING2 . AND . JING . LE . JING3 )  CKS- 
$  CKSS*CKL* (DRS+DRW) / {CKSS*DRW+CKL*DRS) 

ENDIF 

IF(I.EQ.l)  THEN 
RIM-R 
DRl-DR 
CK1M*CKN 
ENDIF 

IF (I. GT. NRS)  THEN 
R-RI- (I-NRS-0 . 5) *DRW 
DR-DRW 

IF ( JING . LT . JING2 . OR . JING . GT . JING3 )  THEN 
CKP-CKSS 
CKN-CKSS 
CKS-CKSS 
CKE-CKSS 
CKW-CKSS 
CKF-CKSS 
CKB-CKSS 
DENP-DENS 
CP-CS 

IFCJING.EQ. (JING2-1))  CKF-2 . 0*CKSS*CKL/ (CKSS+CKL) 
IF(JING.EQ. (JING3+1))  CKB-2 . 0*CKSS*CKL/ (CKSS+CKL) 

ENDIF 

IF ( JING . GE . JING2 . AND . JING . LE . JING3 )  THEN 
CKP-CKL 
CKN-CKL 
CKS-CKL 
CKE-CKL 
CKW-CKL 
CKF-CKL 
CKB-CKL 
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OENP-DBNL 

CP-CL 

IF (I .EQ. (NRS+1) )  CKN-CKSS*CKL* (DRS+DRW) / {CKSS*DRW+CKL*DRS) 
1F(JING.EQ. JING2)  CKB-2 . 0*CKSS*CKL/ (CKSS+CKL) 

1F(JING.EQ. JING3)  CKF-2 . 0*CKSS*CKL/ (CKSS+CKL) 

ENDIF 

ENDIF 


VOLP- ( (R+DR/2) **2 . 0- (R-DR/2) **2 . 0) * (DTHETA/2) *D2 
QHP-QHP+DENP*VOLP*CP* (T (1, J, K) -TSTART) *NGROV 

IF ( I . EQ . NR . AND . JING . NE . JING2 . AND . JING . NE . J1NG3 . AND . 
$  DQDABAR(K) .GT.0.0)  CKS-0.0 


C— CKW*F1 

A-3.0*DENP*CP*DZ**2.0/DT+ (CKW+CKE) *F1 
B— CKE*F1 

D2- (DZ**2 . 0*CKN* (R+DR/2 .0) /R/DR**2 . 0) *F 
D3- (D2**2 . 0*CKS* (R-DR/2 .0) /R/DR**2 . 0) *F 
D4- (DZ**2 . 0*CKB/ (R*DTHETA) **2.0) *F 
D5- (DZ**2 . 0*CKF/ (R*DTHETA) **2 . 0) *F 

IF(J.EQ.l)  D4-0.0 
IF(J.EQ.NTHETA)  D5-0.0 

Dl-3 . 0*DENP*CP*DZ**2 . 0/DT-D2-D3-D4-D5 

IF ( I . EQ . 1 . OR . I . EQ . NR , OR . J . EQ . 1 . OR . J . EQ . NTHEP 1 )  THEN 


D7-(2.0*DZ**2.0*CKS* (R-DR/4.0) /R/DR**2.0) *F 

IF(I.EQ.l.OR.I.EQ.NR)  THEN 
IF(J.EQ.1.0R. J.EQ.NTHEPl)  THEN 
IF(I.EQ.l.AND.J.EQ.l)  D-(D1+D2) *T (1, J, K) +D3*T (I+l, J, K) 

$  +D5*T(I, J+1,K) 

$  + (QWALL ( J, K) *RO*DZ**2 . 0/R/DR) *F 

IF (I. EQ.l. AND. J.EQ.NTHEPl)  D-(D1+D2) *T (I, J, K) +D3*T (I+l , J, K) 

S  +D4*T(I, J-1,K) 

$  + (QWALL ( J,K) *RO*DZ**2 . 0/R/DR) *F 

IFd.EQ.NR.AND.  J.EQ.l)  D- (D1+D3)  *T  (I,  J,  K) +D2*T  (I-l,  J,  K) 

S  +D7*TVOLD(K)  +D5*T(I, J+1,K) -D7*TOLD (I, J,K) 

IFd.EQ.NR.AND.  J.EQ.NTHEPl)  D- (D1+D3)  *T  (I,  J,  K) 

$  +D2*T  (1-1 ,  j,  K)  +D7*TV0LD  (K)  +D4*T  (I ,  J-1 ,  K)  -D7*TOLD  (I ,  J,  K) 

ELSE 

IF  (I.  EQ.l)  D-(D1+D2)*T(I,  J,K)+D3*T(I+1,  J,K)+D4*Td,  J-1,K) 

$  +D5*T (I, J+1 , K) + (QWALL ( J, K) *RO*DZ**2 . 0/R/DR) *F 

IF ( I . EQ. NR)  D- (D1+D3) *T (I , J, K) +D2*T (1-1 , J, K) 

$  +D7*TVOLD{K)+r-4*T(I,  J-1,K)+D5*T(I,  J+l,K)-D7*TOLDd,  J,K) 

ENDIF 
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ELSE 

IF(J.EQ.l)  D=D1*T<I, J,K)+D2*T(I-1, J,K)+D3*T(I+1, J,K) 

$  +D5*T(I, J+1,K) 

IF(J.EQ.NTHEPl)  D-D1*T (I , J, K) +D2*T (I-l , J, K) +D3*T ( I+l , J, K) 
$  +D4*T(I, J-1,K) 

ENDIF 


ELSE 

D-D1*T(I, J,K)+D2*T(I-1, J,K)+D3*T(I+1, J,K)+D4*T(I, J-1,K) 
$  +D5*T(I, J+1,K) 

ENDIF 

IF{D1.LT.0.0)  WRITE(6, 128)1,  J,K 
12S  FORMAT (IX, 'Dl-(-)  IN  ROOP3  I, J, K-* , 3 (2X, 13) ) 
IF(Dl.LT.O.O)  GO  TO  1000 
IF(K.EQ.l)  THEN 
C-0.0 

A-A-CKW*F1 

B-B 

D-D 

ENDIF 

IF(K.EQ.N2)  THEN 

c-c 

A-A-CKE*F1 

B-0.0 

D-D 

ENDIF 

C 

C  TDMA  ALGORITHM 
C 

P(K)— B/(A+C»P(K-1)  ) 

Q(K)-(D-C*Q(K-1)  )  /  (A+C*P(K-1)  ) 


IF(I.EQ.NR)  ^KNRf J,K)-CKS 
IF(K.LE.KBURN)  CKNR ( J, K) -0 . 0 

IF(NDT.EQ.l.AND.I.EQ.NR)  THEN 
QVIN-CKNR ( J, K) * (DTHETA* (RW+DRW/4 .0) *DZ) *NGROV 
$  * (T (NR, J,K)-TVOLD(K))/ (DRW/2.0) 

QV  (K)  -QV  (K)  +(3VIN 

DQVDA(K)-QV(K) / (2 . 0*3 . 1416*RW*DZ) 

ENDIF 

100  CONTINUE 

DO  170  K1-1,NZ 
K-NZ+l-Kl 

TT(I, J,K)-P (K) *TT(I, J,K+1)+Q(K) 

170  CONTINUE 
160  CONTINUE 
180  CONTINUE 


146 


IF{NDT.EQ.l)  GO  TO  182 


AEVAP»NEVAP*DZ* (2 . 0*3 . 1416*RO) 

ACOND- (NZ-NADIA) *DZ* (2 . 0*3 . 1416*RO) 
TCAVG-0 . 0 

DO  110  K-NAD1A+1,NZ 
DO  115  J-1,NTHETA 

TCAVG-TCAVG+TW ( J, K) / (NTHETA* (NZ-NADIA) ) 
115  CONTINUE 
110  CONTINUE 

QWTES^O . 0 
DO  187  K-1,NZ 
QWTES-QWTES+QVTES (K) 

187  CONTINUE 


XXX-AEVAP  *  QEVAP  *  PULSE 

qettl-qettl+aevap*qevap*pulse*dt 

QCTTL-QCTTL+ACOND*EMISS*SlGMA* (TCAVG**4 . 0-TAMB**4 . 0) *DT 
QTES-QTES -QWTES  *DT 


182  DO  193  K-1,NZ 

DO  194  J«1,NTHEP1 
DO  195  I-1,NR 
T(I,  J,K)  -TT{I,  J,K) 
195  CONTINUE 
194  CONTINUE 
193  CONTINUE 


C 

C  &&&&&&&&  COMPUTATION  FOR  THE  SECOND  ONE-THIRD  TIME  STEP  &&&&&&£ 

NADI-123 

DO  280  I-1,NR 
DO  260  K-1,N2 
DO  200  J-1,NTHEP1 


NGORD- ( J-1 ) /NTHEPG+1 
JING-J- <NGORD-l) *NTHEPG 


IF(K.LE.KBURM)  THEN 
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CKL»0. 0000001 
ELSE 

CKL-75.7 

ENDIF 


IF(I.LE.NRS)  THEN 
R«RO- (1-0.5) *DRS 
DR-DRS 
CKP-CKSS 
CKN-CKSS 
CKS-CKSS 
CKE=CKSS 
CKW-CKSS 
CKF-CKSS 
CKB-CKSS 
DENP=DENS 
CP-CS 

IF ( I . EQ . NRS . AND . JING . GE . JING2 . AND . JING . LE . JING3 )  CKS- 

CKSS*CKL* (DRS+DRW) / (CKSS*DRW+CKL*DRS) 

ENDIF 

IF (I. GT, NRS)  THEN 
R-RI- (I-NRS-0 . 5) *DRW 
DR-DRW 

IF (JING . LT . JING2 , OR . JING . GT . JING3 )  THEN 
CKP-CKSS 
CKN-CKSS 
CKS-CKSS 
CKE-CKSS 
CKW-CKSS 
CKF-CKSS 
CKB-CKSS 
DENP-DENS 
CP-CS 

IF(JING.EQ. (JING2-1))  CKF-2 . 0*CKSS»CKL/ (CKSS+CKL) 

IF (JING. EQ. (JING3+1) )  CKB-2 . 0*CKSS*CKL/ (CKSS+CKL) 

ENDIF 


IF ( JING, GE, JING2 . AND . JING. LE . JING3)  THEN 
CKP-CKL 
CKN-CKL 
CKS-CKL 
CKE-CKL 
CKW-CKL 
CKF-CKL 
CKB-CKL 
DENP-DENL 
CP-CL 

IFd.EQ.  (NRS+1)  )  CKN-CKSS^CKL*  (DRS+DRW)  /  (CKSS*DRW+CKL*DRS) 
IF(JING.EQ. JING2)  CKB-2 . 0*CKSS*CKL/ (CKSS+CKL) 

IF(JING.EQ. JING3)  CKF-2 . 0*CKSS*CKL/ (CKSS+CKL) 


ENDIF 

ENDIF 


IF ( 1 . EQ . NR . AND . JING . NE . JING2 . AND . JING . NE . JING3 . AND . 
$  DQDABAR(K) .GT.0.0)  CKS-0.0 


C— CKB*F1 

A-3. 0*DENP*CP* (R*DTHETA) ^*2. 0/DT+ (CKB+CKF) *F1 
B— CKF*F1 

D2-(R*DTHETA**2.0*CKN* (R+DR/2.0) /DR**2.0) *F 
D3-(R*DTHETA**2.0*CKS* (R-DR/2.0) /DR**2.0) *F 
D4- ( (R*DTHETA) **2 . 0*CKW/DZ**2 . 0) *F 
D5«{ (R*DTHETA) **2 . 0*CKE/DZ**2 . 0)  *F 
IF(K.EQ.l)  D4-0.0 
IF(K.EQ.NZ)  D5-0.0 

Dl-3 . 0*DENP*CP* (R^DTHETA) **2 . 0/DT-D2-D3-D4-D5 

IFd.EQ.l.OR.I.EQ.NR.OR.K.EQ.l.OR.K.EQ.NZ)  THEN 

D7-(2.0*R*DTHETA**2.0*CKS» (R-DR/4. 0) /DR**2. 0) *F 

IF(I.EQ.l.OR.I.EQ.NR)  THEN 
IF(K.EQ.l.OR.K.EQ.NZ)  THEN 

IF(I.EQ.l.AND.K.EQ.l)  D-(D1+D2) *T (I, J, K) +D3*T (I+l, J, K) 

$  +D5»T  (I ,  J,  K+1 )  +  (QWALL  ( J,  K)  *RO*R*DTHETA’‘*2 . 0  /DR)  *F 

IF(I.EQ.l.AND.K.EQ.NZ)  D“(D1+D2) *T (I, J, K) +D3*T (I+l , J, K) 

$  +D4*T (I, J, K-1 )  +  (QWALL ( J,  K) *RO*R*DTHETA**2 . 0/DR) *F 

IFd.EQ.NR.AND.K.EQ.l)  D-(D1+D3-D7)  *T  (I,  J,  K) +D2*T  (1-1 ,  J,  K) 

S  +D7*TVOLD(K)+D5*T(I, J,K+1) 

IF(I.EQ.NR.AND.K.EQ.NZ)  D“(D1+D3-D7) *T (I , J, K) +D2»T (1-1 , J, K) 

$  +D7*TVOLD(K)+D4*T(I, J,K-1) 

ELSE 

IFd.EQ.l)  D-(D1+D2)*T(I,  J,K)+D3*T(I+1,  J,K)+D4*T(I,  J,K-1) 

$  +D5*T (I, J, K+1 )  +  (QWALL ( J, K) *RO*R*DTHETA*  *2 . 0/DR) *F 

IF (I .EQ.NR)  D- (D1+D3-D7) *T (I, J,K) +D2*T (1-1, J, K) 

$  +D7*TVOLD(K)  +D4*T (I, J, K-1 ) +D5*T (I, J, K+1) 

ENDIF 

ELSE 

IF(K.EQ.l)  D-D1*T(I,  J,K)+D2*Td-l,  J,K)+D3*T(I+1,  J,K) 

$  +D5*T(I, J,K+1) 

IF(K.EQ.NZ)  D-D1*T(I,  J,K)+D2*Td-l,  J,K)+D3*T(I+1,  J,K) 

$  +D4*T(I,  J,K-1) 

ENDIF 


ELSE 

D-D1*T(I,  J,K)+D2*Td-l,  J,K)+D3*T(I  +  1,  J,K)+D4*T(I,  J,K-1) 
$  +D5*T(I, J,K+1) 

ENDIF 
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224  IF(Dl.LT.O.O)  WRITE (6, 228) 1,  J,  K 
228  FORMAT (IX, 'Dl-C-)  IN  ROOP2  I, J, K- ' , 3 (2X, I3) ) 
IF(Dl.LT.O.O)  GO  TO  1000 
IF(J.EQ.l)  THEN 
A-A-CKB*F1 
B-B 
D-D 
END  IF 

IF(J.EQ.NTHEPl)  THEN 
C=C 

A-A-CKF*F1 

D=D 

ENDIF 

C 

CTDMA  ALGORITHM 
C 

IF(J.EQ.l)  THEN 
P(J)— B/A 
Q(J)-D/A 
ELSE 

P(J)=-B/ (A+C*P(J-1)) 

Q(J)-(D-C*Q(J-1) ) / (A+C*P{J-1) ) 

ENDIF 

200  CONTINUE 

DO  270  J1-1,NTHEP1 
J-NTHEPl+l-Jl 

IF(J.EQ.NTHEPl)  TT (I, J, K) -Q (J) 

IF ( J . NE . NTHEPl )  TT ( I ,  J,  K) -P ( J) *TT ( I , J+1 , K) +Q ( J) 
270  CONTINUE 
260  CONTINUE 
280  CONTINUE 

DO  293  K-1,NZ 
DO  294  J-1, NTHEPl 
DO  295  I-1,NR 
T(I,  J,K)-TT(I,  J,K) 

295  CONTINUE 
294  CONTINUE 
293  CONTINUE 


IF{IPCM.EQ,2)  THEN 
CALL  TES 
ENDIF 

C 

C  &&&&&&&£  COMPUTATION  FOR  THE  THIRD  ONE-THIRD  TIME  STEP  &&&&&&&& 
NADI-  4 

NQVO-0 

350  DO  405  K-1,NZ 
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n  o 


QVOLD (K) -QV (K) +QVTES (K) 
TVOLD(K)-TV(K) 

405  CONTINUE 


NQVO-NQVO+1 

Cl  1F(NQV0.GT.  1)  GO  TO  450 


DO  448  K-1,NZ 

IF{NQV0.GT.1)  DQVDA(K)-DQVDA(K)+FRELAX* (DQVDAN(K) -DQVDA(K) ) 

IF  <NQV0 . GT . 1 )  DQDATE {K) "DODATE (K) +FRELAX* (DQDATN (K) -DQDATE (K) ) 
448  CONTINUE 


GUESS  TV(KBURN+1)  FOR  THE  NEXT  TIME  STEP 

TVEEl-0.0 
TVEE2-0 . 0 
NCALLV-0 

400  NCALLV-NCALLV+l 

IF(NCALLV.GT.  200)  THEN 
WRITE(6,») 'NCALLV  >  200' 

GO  TO  1000 
ENDIF 


I F ( NCALLV . EQ , 1 )  TV ( KBURN+ 1 ) »TV0 ( KBURN+ 1 ) +DTV1 
I F ( NCALLV . EQ . 2 )  TV ( KBURN+ 1 ) -TVO ( KBURN+ 1 ) -DTVl 
IF (NCALLV. GT. 2)  THEN 
IF( (TVEE1*TVEE2) .EQ.O.O)  THEN 
WRITE (6, *) 'TVEE1*TVEE2»0.0' 

GO  TO  1000 
ELSE 

TV(KBURN+1)-(TVEE1+TVEE2) /2.0 
ENDIF 
ENDIF 


C  IF( (NDT/XWRITE) .EQ. (NDT/NWRITE) )  THEN 

C  WRITE ( 6, * )  ' NCALLV- ' , NCALLV, • TVO- * , TV (KBURN+1 ) 
C  ENDIF 


Cl  CALL  BOWMAN 
CALL  CHI 


DO  415  N-  1,NZ 

Cl  TV(N)-TV(KBURN+1) 
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TV(N)-TVREAL(N) 

415  CONTINUE 


DO  380  K-1,NZ 
DO  360  J-1,NTHEP1 
NGORD- (J-1) /NTHEPG+1 
JING-J- (NGORD-1) ‘NTHEPG 
DO  300  I-1,NR 


IF(K.LE.KBURN)  THEN 
CKL-0. 0000001 
ELSE 

CKL-75.7 

ENDIF 


IF(I.LE.NRS)  THEN 
R-RO- (I-O . 5) *DRS 
DR-DRS 
CKP-CKSS 
CKN-CKSS 
CKS-CKSS 
CKE-CKSS 
CKW-CKSS 
CKF-CKSS 
CKB-CKSS 
DENP-DENS 
CP-CS 

IFd.EQ.NRS.AND.  JING.GE.  JING2.AND.  JING.LE.  JING3)  CKS- 
$  CKSS*CKL* (DRS+DRW) / (CKSS*DRW+CKL*DRS) 

ENDIF 


IF(I.GT.NRS)  THEN 
R-RI- (I-NRS-0 .5) *DRW 
DR-DRW 

IF ( JING . LT . JING2 . OR . JING . GT . JING3 )  THEN 
CKP-CKSS 
CKN-CKSS 
CKS-CKSS 
CKE-CKSS 
CKW-CKSS 
CKF-CKSS 
CKB-CKSS 
DENP-DENS 
CP-CS 

IF(JING.EQ. (JING2-1))  CKF-2 . 0*CKSS*CKL/ (CKSS+CKL) 
IF (JING. EQ. {JING3+1))  CKB-2 . 0*CKSS*CKL/ (CKSS+CKL) 
ENDIF 


152 


IF ( JING . GE . J1NG2 . AND . JING . LE . JING3 )  THEN 
CKP-CKL 
CKN-CKL 
CKS-CKL 
CKE-CKL 
CKW-CKL 
CKF-CKL 
CKB-CKL 
DENP-DENL 
CP-CL 

IFd.EQ.  (NRS+D)  CKN=CKSS*CKL*  (DRS+DRW)  /  (CKSS*DRW+CKL*DRS) 
IF ( JING. EQ. JING2)  CKB=2 . 0*CKSS*CKL/ (CKSS+CKL) 

IF(JING.EQ.  JING3)  CKF-2 . 0*CKSS*CKI,/  (CKSS+CKL) 

ENDIF 
END  IF 


IF ( I . EQ . NR . AND . JING . NE . JING2 . AND . JING . NE . JING3 . AND . 
S  DQDABAR(K) .GT.0.0)  CKS=0.0 


C— CKN*  (R+0 . 5*DR)  *F1 

A-3 . 0*DENP*CP*R*DR**2 . 0/DT+ (CKN* (R+0 . 5*DR) 

$+CKS* (R-0 . 5*DR) ) 'FI 
B— CKS*  (R-0 . 5*DR)  *F1 
D2- (DR**2 . 0»CKB/R/DTHETA**2 . 0) *F 
D3- (DR**2 . 0*CKF/R/DTHETA*»2 .0) *F 
D4= (R*DR**2. 0*CKW/DZ**2. 0) *F 
D5- (R*DR**2 . 0*CKE/DZ**2 . 0) *F 

IF(J.EQ.l)  D2-0.0 
IF(J.EQ.NTHETA)  D3»0.0 
IF(K.EQ.l)  D4-0.0 
IF(K.EQ.NZ)  D5-0.0 

Dl-3 . 0*DENP*CP»R*DR**2 . 0/DT-D2-D3-D4-D5 
IF ( J . EQ . 1 . OR . J . EQ . NTHEPl . OR . K . EQ . 1 . OR . K . EQ . NZ )  THEN 
IF(J.EQ.l.OR.J.EQ.NTHEPl)  THEN 
IF(K.EQ.l.OR.K.EQ.NZ)  THEN 
IF(J.EQ.l.AND.K.EQ.l)  D-D1*T (I , J, K) 

$  +D3*T(I, J+1,K)+D5*T(I, J,K+1) 

IF (J.EQ. 1 .AND.K.EQ.NZ)  D-D1*T (I, J, K) 

$  +D3*T(I, J+1,K)+D4*T(I, J,K-1) 

IF (J.EQ. NTHEPl .AND. K.EQ. 1)  D-Dl *T ( I , J, K) +D2*T ( I , J-1 , K) 

S  +D5*T(I, J,K+1) 

IF (J.EQ. NTHEPl .AND.K.EQ.NZ)  D=D1 *T ( I , J, K) +D2*T ( I , J-1 , K) 

$  +D4*T(I, J,K-1) 

ELSE 

IF(J.EQ.l)  D-D1*T(I, J,K)  +D3*T ( I , J+1 , K) 

S  +D4*T(I,J,K-1)+D5*T(I, J,K+1) 

IF (J.EQ. NTHEPl)  D=D1*T (I , J, K) +D2*T ( I , J-1 , K) 
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+D4*T{1, J,K-1)+D5*T(I, J,K+1) 


ENDIF 

ELSE 

IF(K.EQ.l)  D=D1*T(1, J,K)+D2*T(I, J-1,K)+D3*T(I, J+1,K) 

$  +D5*T(1, J,K+1) 

IF(K.EQ.NZ)  D=D1*T(I, J,K)+D2*T(I, J-1,K)+D3*T(I, J+1,K) 

S  +D4*T<I, J,K-1) 

ENDIF 

ELSE 

D-D1*T(I, J,K)+D2*T(I, J-1,K)+D3*T(I,  J+1,K)+D4*T(I,  J,K-1) 
$  +D5*T(I, J,K+1) 

ENDIF 

IF(Dl.LT.O.O)  WRITE(6, 328)1,  J,K 
328  FORMAT (IX, 'Dl-(-)  IN  ROOPl  I , J, K=' , 3 (2X, 13) ) 
IF(Dl.LT.O.O)  GO  TO  1000 

IF(I.EQ.l)  THEN 
C=0.0 

A-A-CKN* (R+DR/2 . 0) *F1 
B=B 

D=D+ (QWALL ( J, K) *RO*DR) *F1 
ENDIF 

IF(I.EQ.NR)  THEN 
C-C 

A=A+CKS*R*F1 

B*0.0 

D-D+2 . 0*CKS*  <R-DR/4 . 0)  *.TV  (K)  *F1 
ENDIF 
C 

C  TDMA  ALGORITHM 
C 

P(I)— B/(A+C*P(I-1)  ) 

Q(I)=<D-C*Q(I-1) ) / (A+C*P(I-1) ) 


300  CONTINUE 

DO  370  I1-1,NR 
I-NR+l-Il 

TT(I,  J,K)-P(I)  *TT(I+1,  J,K)+Q(I) 
370  CONTINUE 
360  CONTINUE 
380  CONTINUE 


IF(IPCM.EQ.2)  THEN 
CALL  TES 
ENDIF 
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o  o  o  o  o  o  o 


SEARCHING  FOR  TV(KBURN+1)  FOR  NEXT  TIME  STEP 

QIN-0.0 
SUMQVI=0 . 0 
SUMQVO-0 . 0 

DO  410  K»1,NZ 
QV(K)*0.0 
DO  420  J-1,NTHETA 
NGORD- ( J-1) /NTHEPG+1 
JING-J- (NGORD-1) *NTHERG 

DO  430  I-1,NR 

IF (I -EQ.NR)  QVIN=CKNR(J,K) * (DTHETA* (RW+DRW/4 . 0) *DZ) *NGROV 
$  *(TT(NR,J,K)-TV(K>)/ (DRW/2.0) 

430  CONTINUE 

QV(K)-QV(K)+QVIN 
420  CONTINUE 

IF (QV (K) , GT . 0 . 0 )  QIN-QIN+QV (K) 

IF< (QV{K)+QVTES (K) ) .GT.0.0)  SUMQV1-SUMQVI+ (QV (K) +QVTES (K) ) 
IF  (  (QV  (K)  +QVTES  (K)  )  .  LE.  0 . 0)  SUMQVO-SUMQVO-  (QV  (K)  +QVTES  (K)  ) 


410  CONTINUE 

SUMQV-SUMQVI -SUMQVO 

IF( (NDT/XWRITE) .EQ. (NDT/NWRITE) )  THEN 

WRITE(6,*)  ' SUMQV= • , SUMQV, 'QVI=' , SUMQVI, •QVO=' , SUMQVO 

ENDIF 

IF ( (ABS (SUMQV) /ABS ( S UMQV I + SUMQVO) ) .GT.EPSV)  THEN 
IF ( SUMQV. LE. 0.0)  TVEEl-TV (KBURN+1 ) 

IF (SUMQV. GT. 0.0)  TVEE2 -TV (KBURN+1 ) 

GO  TO  400 
ELSE 

GO  TO  440 
ENDIF 


440  IF ( (NDT/XWRITE) .EQ. (NDT/NWRITE) )  THEN 
C  WRI TE ( 6 , * ) ' NCALLV- • , NCALLV 

ENDIF 


Cl  CALL  BOWMAN 

Cl  CALL  CHI 


PVMIN-PV (KBURN+1 ) 
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DO  465  N-KBURN+2,NZ 
IF (PV (N) . LT . PVMIN)  THEN 
PVMIN-PV(N) 

KPVMIN-N 
END  IF 

465  CONTINUE 


CALL  LIQUID 


IF( (NDT/XWRITE) .EQ. (NDT/NWRITE) )  THEN 
WRITE(6,*)'  Z  PV  PL  DPV 

DO  460  N-KBURN+1,NZ 
Z=(ZL/NZ) * (N-0.5) 

DPV(N)-PV(N) -PV(1) 

DPL(N)-PL(N) -PL(1) 

WRITE (6, 462)  Z, PV(N) , PL (N) , DPV (N) , DPL (N) 

462  FORMAT<2X,F8.4,4 (2X,F8.1) ) 

460  CONTINUE 
ENDIF 


DPVMAX-PV (KBURN+1 ) -PVMIN 
DPLMAX-PL(NADIA+1) -PL(KBURN+1) 

PCMAX-2. 0*TENLIQ/BC 

IF ( (NDT/NWRITE) .EQ. (NDT/XWRITE))  WRITE(6,*) ' DPVMAX- ' , DPVMAX 
$  , 'DPLMAX-',DPLMAX, 'PCMAX-',PCMAX 


DO  442  K-1,NZ 
QVNEW  (K)  -QV  (K)  +QVTES  (K) 

TVNEW(K)-TV(K) 

DQVDAN(K)-QV(K) / (2 . 0*3 . 1416*RW*DZ) 
IF(lPCM.EQ.l)  GO  TO  442 

DQDATN (K) -QVTES (K) / (NTESTU*2 . 0*3 . 1416*RTES*DZ) 
442  CONTINUE 


IF (METHOD. EQ.l)  GO  TO  450 


C  DO  445  K-1,NZ 

C  WRITE(6,906)K,QV(K) , QVTES (K) , QV (K) +QVTES (K) 

C  IF(K.GT.l .AND. ( (QV (K) +QVTES (K) ) * (QV (K-1 ) +QVTES (K-1 ) ) ) . LT.O. 0) 

C  $PLM-K*DZ 
C445  CONTINUE 

C  IF ( (NDT/XWRITE) .EQ. (NDT/NWRITE) )  THEN 

Cl  WRITE(6, *) 'PLM-' ,PLM 


DPL' 
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Cl  WRITE (6,*)'  ' 

Cl  WRITE (6,*)'  • 

Cl  WRITE(6,’*)  'NQV0-',NQV0 

Cl  WRITE (6, *) 'QVO-' ,QVO (1) , * QVOOLD- * , QVOLD ( 1 ) , 'QVONEW-' , QVNEW ( 1 ) 
Cl  WRITE (6, *) 'TVO-' ,TV0 (1) , 'TVOOLD-' , TVOLD (1 ) , ' TVONEW- * , TVNEW ( 1 ) 
Cl  WRITE (6,*)'  ' 

Cl  WRITE (6,*)'  ' 

C  ENDIF 


Cl  IF(ABS ( (TVNEW(l) -TVOLD(l) ) / (TVNEW(l) -TVO (1) ) ) .LE.EPSQVK)  THEN 

Cl  WRITE (6,*)'  ' 

Cl  WRITE (6,*) 'TVO  CONVERGE  NOW!' 

Cl  WRITE (6,*)'  ' 

Cl  ENDIF 


NWOO-0 

DO  452  K-1,NZ 

IF  (ABS  (  (TVNEW  (K)  -TVOLD  (K)  )  /  (TVNEW  (K)  -TVO  (K)  )  )  .  GT .  EPSQVK) 
$  NWOO-NWOO+1 

IF(NWOO.GE.30)  GO  TO  350 

452  CONTINUE 


450  DO  455  K=1,NZ 

DQVDA (K) -DQVDAN (K) 
DQDATE (K) -DQDATN (K) 
455  CONTINUE 


NADI-6 

IF(IPCM.EQ.2)  THEN 
CALL  TES 
ENDIF 

THPTTL-0 . 0 
TVTTL-0 . 0 
DO  393  K-1,NZ 
DO  394  J-1,NTHEP1 
NGORD- ( J-1 ) /NTHEPG+1 
JING-J- (NGORD-1 ) *NTHEPG 
DO  395  I-1,NR 

IF ( I . GT . NRS . AND . JING . GE . JING2 . AND . JING . LE . JING3 . AND . K . LE . KBURN) 

$  THEN 

T(I,  J,K)-TV(K) 

ELSE 

T{I, J,K)-TT(I, J,K) 
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ENDIF 

TOLDd,  J,K)=T(I,  J,K) 

THPTTL-THPTTL+T (I , J, K) 

395  CONTINUE 

TW(J,K)-T(1, J,K)+(DRl/2. 0) *QWALL(J,K) *R0/CK1M/ (RlM+DRl/4 . 0) 
394  CONTINUE 

TVTTL-TVTTL+TV (K) 

DTV (K) =TV (K) -TVOLD (K) 

393  CONTINUE 


THPAVG (NDT) -THPTTL/NZ/NTHEPl/NR 
TVAVG (NDT) -TVTTL/NZ 
CTIME (NDT) -TIME 

IF ( (NDT/XWRITE) . EQ . (NDT/NWRITE) )  THEN 

WRITE(6,910) 

910  FORMAT (IX, //////) 

WRITE (6, 920) TIME 

920  FORM?*T(28X,  '  TIME- ' ,  2X,  F7 . 2, 2X,  'SECONDS',/) 


WRITE (6,*)'  ' 

WRITE (6,*)'  ' 

WRITE (6,*)  '  ' 

WRITE(6,940)  (TW(2,K) ,K=1,NZ,  4) 

940  FORMATdX,  '  WALL' ,  10  ( IX,  F6 . 1)  ) 

DO  900  I-1,NR 

WRITE (6, 950) I, (T(I,2,K),K=1,NZ,  4) 

950  FORMATdX, 'I-',  13, 10  (IX,  F6.1)) 

900  CONTINUE 

WRITE(6,960)  (TV(K) ,K-1,NZ,  4) 

960  FORMATdX,  'VAPOR' ,  10  (1X,F6. 1) ///) 

WRITE (6, 970)  THPAVG (NDT) 

970  FORMAT ( 2 X, 'AVERAGE  HEAT  PIPE  TEMPERATURE- ', 2X, F9 . 4 ) 
WRITE (6, 980)  TVAVG (NDT) 

980  FORMAT (2X, 'AVERAGE  VAPOR  TEMPERATURE- ', 2X, F9 . 4/// ) 

ENDIF 

C 

C 

C 


IF ( (NDT/XWRITE) . EQ . (NDT/NWRITE) )  THEN 
IF(IPCM.EQ.2)  THEN 

WRITE (6,*)'  TES  TEMPERATURE  ' 

WRITE (6,*)'  ' 
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DO  990  I-1,NRTES,2 

WRITE(6,930)  I, (TE<1,K) ,K=1,NZ,  4) 

930  FORMAT (IX, ' IE»' , 13, 10 (1X,F6. 1) ) 

990  CONTINUE 
ENDIF 
ENDIF 

IF ( (NDT/NWRITE) . EQ . (NDT/XWRITE) )  THEN 
PCTQWTES=(QWTES/QIN) *100.0 
WRI TE ( 6 , * ) • SUMQVI  =  * , SUMQVI , ’ SUMQVO- ' , SUMQVO 
WRITE (6,*)  'QIN-',QIN, ' QWTES« * , QWTES , ' %QWTES- ' , PCTQWTES, '%' 
WRITE ( 6, * ) ' QEVAP= ' ,  XXX 
WRITE (6,*)'  ' 

WRITE (6,*)' 


C  ENERGY  BALENCE  CHECK 

IF(NDT.GT.1.AND. ( (NDT/NWRITE> .EQ. (NDT/XWRITE) ) )  THEN 
ERRORl= (QETTL- (QCTTL+QHP+QPCM) ) /QETTL*100 . 0 
ERROR2” (QETTL- (QCTTL+QHP+QTES) ) /QETTL*100 . 0 
IF(IPCM.EQ.2)  ERROR3- (QTES-QPCM) /QTES*100.0 
IF ( IPCM . EQ . 2 )  PCTMELT-QMELT/QPCMTTL* 100.0 

WRITE (6, 902)  QETTL, QCTTL,QHP,QPCM,QTES 

902  FORMATdX,  'QETTL- F9 . 1, 2X,  *  QCTTL-' ,  F9 . 1 , 2X,  '  QHP-’ ,  F8 . 1 , 2X 

$  , 'QPCM-' ,F8.1,2X, 'QTES-' ,F8.1/) 

WRITE (6,903)  ERRORl , ERROR2 , ERROR3 , PCTMELT 

903  FORMATdX, 'ENERGY  BALENCE  ERROR- ',2X,F6. 2,  '  %',2X,F6.2,'  %' 

$  ,2X,F6.2,'  %',2X, '%MELT=',F6.2,  '  %'///) 

ENDIF 


C  DO  905  K-1,NZ 

C  WRITE (6, 906) K, QV (K) , QVTES (K) , QV (K) +QVTES (K) 

C  906  FORMAT (IX, 13, 3X, 'QV-' ,F10. 5, 5X, ' QVTES-' , FIO . 5, 5X, 'QVNET-' , FIO . 5) 
C  905  CONTINUE 
C  WRITE(6,*) 'PLM«',PLM 

C  WRITE (6,*) 'NQV0-',NQV0, ' FRELAX- ' , FRELAX 

DO  915  N-1,NZ 

WRITE(6,916)  N,  PV(N) ,DPV(N) ,TVREAL(N) ,DENV(N) 

916  FORMATdX,  'N-',I3,2X,  '  PV  (N) -' ,  IX,  F8 . 1, 2X,  'DPV  (N) -' ,  IX,  F8 . 1 , 

$  2X, 'TV(N)-',1X,F7.2,2X, ' DENV (N) -' , IX, F8 . 6) 

915  CONTINUE 


ENDIF 

IF  (NDT.EQ.NDTEND.OR.NQVO.GT. 600)  THEN 

DO  982  K-1,NZ 
WRITE (6, 983)  K,DQVDA(K) 

C983  FORMAT (3X, 'K-',I3,2X, 'DQVDA-' , F15. 3) 
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C982  CONTINUE 


GO  TO  999 
ENDIF 


60  CONTINUE 


999  IF(IPCM.EQ.l)  WRITE{6,*)*C  WITHOUT  TES* 

IF(IPCM.EQ.2)  WRITE(6,*)'C  WITH  TES' 

IF (METHOD. EQ.l)  WRITE(6,*)'C  EXPLICIT  WITH  VAPOR' 
IF(METHOD.EQ.2)  WRITE(6,*)'C  IMPLICIT  WITH  VAPOR' 

WRITE(6, *) 'C , ' NTESTU« ' , NTESTU, ' RTES-' , RTES, 'DT«' ,DT, 'F=' ,F 
WRITE(6,*)'C  TIME  -THP-  -TV-' 

DO  995  NDT«NWRITE,NDTEND,NWRITE 

WRITE (6, 996) CTIME (NDT) , THPAVG (NDT) , TVAVG (NDT) 

996  F0RMAT(1X,F7.2,2  (2X,F8.2) ) 

995  CONTINUE 

1000  STOP 
END 

C 

C 

C 

C 

^  it  it  it  ii  Hit  it  it  •ti itit  itit  it  •ik  itiri/t 

C  *  SUBROUTINE  BOWMAN  * 

^  it -tt  it  it  iHr  it  ir  it  it  ink  f  ir  it  it  ir  ir  ^t  it 

c 

c 

SUBROUTINE  BOWMAN 

C  ONE  DIMENSIONAL  COMPRESSIBLE  VAPOR  FLOWFOR  HEAT  PIPE 
REAL  *8  TV(IOO) 

COMMON  /HPIPE/  TV,DQDATE{100) ,PL,DZ,NZ,NDT,QVTES (100) , 
SNTESTU, TIME, KBURN, NWRITE, XWRITE, NGROV, RW, NADIA, PV (100 ) 
COMMON  /VCORE/  DQDA(IOO) ,DQDABAR(100) , IPCM,PVTTL(100) 

$  ,TVREAL(100) ,PLM,DTES,DVAPOR,DWICK,DHYDRO,XMMAX,KPLM 
$  ,NEVAP,XTI,DENV(100) 

IF(IPCM.EQ.l)  THEN 
DO  10  N-1,NZ 
QVTES (N)«0.0 
DQDATE(N)-0.0 
10  CONTINUE 
ENDIF 


C  GUESS  THE  VAPOR  TEMPERATURE  AT  THE  EVAPORATOR  END  TO 
T0-TV(KBURN+1) 
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c 

OMEGA-0.25 
R=361 . 5 

C  WRITE (*,*) 'RCONST=',R 

lambda-4636437.26-180. 817*T0 
GAMMA-1 . 67 

GAMMl -GAMMA/ (GAMMA-1 .0) 

INITIALIZE  THE  VARIABLES  AT  THE  UPSTREAM  PIPE  END 

W-0.0 
Z-0.0 

RHO-6.335E08*10.0** (-5567. O/TO) /TO* *1.5 
PO-RHO*R*TO 
C 

C  USE  INCOMPRESSIBLE  MODEL  TO  GET  THE  STARTING  VALUES 

C 


RHOV— DQDABAR  (KBURN+1 )  /LAMBDA 
RHOVTE-RHOV 


C2-GAMMA*P0/RHO 

XM22-(4.0* (RHOV*DWICK+NTESTU*RHOVTE*DTES) *DZ/RHO/DVAPOR**2 . 0)  **2 
S  /C2 

W2-W-RHOV*DZ*DWICK*3. 1416-NTESTU*RH0VTE*DZ*DTES*3. 1416 
Z2-Z+DZ/2.0 

P2-P0/ (1.0+ (GAMMA-1. 0)/2.0*XM22)**GAMMl 
P02-P0 
RH02-RH0 
TV (KBURN+1) -TO 
TVREAL (KBURN+1 ) -TO 
PV (KBURN+1) -PO 
PVTTL (KBURN+1) -PO 
DENV (KBURN+1 ) -RHO 
C 

C  MARCHING  DOWN  THE  PIPE  FINDING  THE  FLOW  PROPERITIES  ALONG  THE  WAY 

C 

Z-Z2 

C 

c 

c 

CC  WRITE (6,*)'  Z/PL  P2/P0  T 

CC  $  W2  XM22' 


XMMAX-0 . 0 

DO  5  N-KBURN+2,NZ 
XM2-XM22 

IF(XM2.GT.0.99)  THEN 
WRITE(6,70) Z,XM2 
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70  FORMAT (3X, ' Z- ' , F6 . 3, 3X, ’XM2-',F6.3) 
ENDIF 

P-P2 

P01«P02 

W-W2 

Z=Z+DZ 

C 

C  PREDICTOR  AND  CORRECTOR  STEPS 
C 

C  CALCULATE  THE  INFLUENCE  COEFFICIENTS 
C 

XM2BAR-(XM2+XM22) /2.0 
CONl-1 . 0+ (GAMMA-1 . 0) /2 . 0*XM2BAR 
CON2=1.0-XM2BAR 
CON3=l . 0+GAMMA*XM2BAR 
CON4 -GAMMA*XM2 BAR 
FW-2 . 0*CON3*CON1/CON2 
FWP-CON4 
C 

C  FIND  THE  MASS  ADDED  THROUGH  THE  PIPE  WALL 
C 

RHOV— DQDABAR(N)  /  LAMBDA 
RHOVTE-RHOV 


DELW— RHOV*DZ*DWICK*  3.141 6-NTESTU*RHOVTE*DZ*DTES  *3.1416 
W2-W+DELW 
IF(W2.LE.0.0)  THEN 
XM22=0. 0 
P02-P01 
T-TO 

GO  TO  30 
ENDIF 
C 

C  FIND  THE  FRICTION  FACTOR  AND  THE  FRICTION  INFLUENCE  COEFFICIENT 
C 

T-TO /CONI 

RMU-6. 083E-09*T+1 .2606E-05 
WBAR-(W+W2) /2.0 

IF(WBAR.LE.O.O)  THEN 
XM22-0.0 
P02-P01 
T-TO 

GO  TO  30 
ENDIF 

REY- (4 . 0*WBAR*DHYDRO) / (3 . 1416*DVAPOR**2 . 0*RMU) 

RER— DELW/  (3.1416*RMU*DZ) 

ZBAR- ( (Z-XTI) /OMEGA) **2 
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IF( (NDT/NWRITE) .EQ. (NDT/XWRITE) )  THEN 
C  WRITE (6,*)  ’REY-',REY 

ENDIF 

IF(REY.GE. 12000)  WRITE(6,*) 'REY>12000' 

FL-16/REY* (1. 2337-0. 2337*EXP(-0.036*ABS(RER) ) ) *EXP (1 . 20*XM2) 

C  FL-24/REY 

BETA-ABS (RER/REY) 

FSTAR-0.046/REY**0.2 

FT-FSTAR* (1+55*REY**0.1*EXP(1.2*XM2) *BETA**0. 9* (PL/DHYDRO) **0.1) 
IF(2.LE.XTI)  F=FL 

IF(Z.GT.XTI)  F=FT-(FT-FL) *EXP (-0 . 412*ZBAR) 


FF=CON4  *  CONI / CON2  *  4 . 0  *  F / DHYDkO 
FFP-CON4/2 . 0*4 . 0*F/DHYDRO 
C 

IF(W2.LE.0.0)  WRITE<6,*)  •W2<0.0' 
XM22-EXP (ALOG (XM2) +FF*D2+FW*ALOG (W2/W) ) 
IF(XMMAX.LT. (XM22**0.5) )  XMMAX-XM22**0 . 5 
P02“EXP (ALOG(POl) -FFP*DZ-FWP*ALOG(W2/W) ) 
30  P2-P02/ <1.0+ (GAMMA-1 .0) /2*XM22) **GAMMl 

TX-T 

DO  80  NX»1,100 

T2— 5567  .  O/ALOGIO  (P2*TX**0 . 5/2 . 29E11) 
IF(ABS(T2-TX) .LT.1.0)  GO  TO  90 
TX-T2 

IF(NX.EQ.IOO)  WRITE{6,*)  'NX=100' 

80  CONTINUE 

90  TVREAL(N)-T2 
PV(N)-P2 
PVTTL(N)-P02 

DENV  (N)  -PV  (N)  /  (R*TVREAL  (N)  ) 

C 

CC  WRITE (6, 20) Z/PL, P2/P0, T, W2, XM22 

CC  20  FORMAT(1X,5(3X,F15.10) ) 

5  CONTINUE 


DO  50  N-1,KBURN 
TVREAL (N) -TVREAL (KBURN+1 ) 
PV(N)-PV(KBURN+1) 

PVTTL (N) -PVTTL (KBURN+1 ) 

50  CONTINUE 
40  RETURN 
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END 


C 

C 

C 

Q  ******************** 

C  *  SUBROUTINE  CHI  * 

Q  ★*★*★**★★★*★★****★★* 

C 

C 

SUBROUTINE  CHI 

C  ONE  DIMENSIONAL  COMPRESSIBLE  VAPOR  FLOW  FOR  HEAT  PIPE 

REAL  *8  TV(IOO) ,W(100) 

REAL  LAMBDA 

COMMON  /HPIPE/  TV,DQDATE{100) ,PL,DZ,NZ,NDT,QVTES (100)  , 
SNTESTU, TIME, KBURN, NWRITE, XWRITE, NGROV, RW, NADIA, PV (100 ) 
COMMON  /VCORE/  DQDA(IOO) ,DQDABAR(100) , IPCM,PVTTL(100) 

$  ,TVREAL(100) ,PLM,DTES,DVAPOR,DWICK,DHYDRO,XMMAX,KPLM 
$  ,NEVAP,XTI,DENV(100) 

DIMENSION  PP(IOO) 

IF(IPCM.EQ.l)  THEN 
DO  10  N=*1,NZ 
QVTES (N) -0 . 0 
DQDATE (N) -0 . 0 
10  CONTINUE 
ENDIF 

AV-3.1416*DVAPOR**2.0/4.0 


C  GUESS  THE  VAPOR  TEMPERATURE  AT  THE  EVAPORATOR  END  TO 
T0=TV(KBURN+1) 

RHO«6.335E8*10.0** (-5567/T0) /T0**1.5 
P0-RHO*R*T0 

R-361 . 5 

LAMBDA-4636437.26-180.817*T0 
GAMMA-1 . 67 

GAMMl -GAMMA/ (GAMMA-1 .0) 

C2-GAMMA*P0/RHO 


W(KBURN+l)-0.0 
TVREAL (KBURN+1 ) -TO 
PV(KBURN+1)-P0 
PP (KBURN+1 )-P0 
DENV (KBURN+1 ) -RHO 
XMMAX-0 . 0 
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DO  5  N-KBURN+2,N2+1 


RMU-6 . 083E-9*TVREAL (N-1) +1 . 2606E-5 

RHO-6.335E8*10.0** (-5567/TVREAL (N-1) ) /TVREAL (N-1) **1 . 5 


C  FIND  THE  MASS  ADDED  THROUGH  THE  PIPE  WALL 
C 

RHOV- -DQDABAR ( N - 1 ) / LAMBDA 
RHOVTE-RHOV 


DELW— RHOV*D2*DWICK*3. 1416-NTESTU*RHOVTE*DZ*DTES*3. 1416 
W(N)-W(N-1)+DELW 


IF(W(N) .LE.0.0)  THEN 
PP(N)-PP(N-1) 
PV(N)-PV(N-1) 

TVREAL (N) -TVREAL (N-1 ) 
GO  TO  5 
ENDIF 


XM2-W(N) **2.0/ (AV*RHO)**2.0/C2 
1F(XMMAX.LT. (XM2**0.5) )  XMMAX-XM2  **0.5 
REY-4.0*DHYDRO*W(N) / (3 . 1416*DVAPOR**2 . 0*RMU) 

IF( (NDT/NWRITE) .EQ. (NDT/XWRITE) )  THEN 
C  WRITE(6,*)  'REY-',REY 

ENDIF 

A-32 . 0*RMU/ (DHYDRO**2 . 0*AV*RHO*LAMBDA) 

B-0 . 076*RMU/ (DHYDRO**2 . 0*AV*RHO* LAMBDA) 

C-1 . 0+ (GAMMA-1 ) /2*XM2 

IF (REY . LE . 2300 . 0 . AND . XM2 . LE . 0 . 04 )  FV-A 
IF(REY. LE. 2300.0. AND. XM2.GT. 0.04)  FV-A*C** (-0.5) 

IF (REY. GT. 2300.0. AND. XM2.LE. 0.04)  FV=B*REY**0 . 75 

IF (REY. GT. 2300.0. AND. XM2.GT. 0.04)  FV-B*REY**0 . 75*C** (-0.75) 

IF (REY. LE. 2300.0)  BETA-1.33 
IF (REY. GT. 2300.0)  BETA-1.0 

PP (N) -PP (N-1 ) -FV*W (N) *LAMBDA*DZ 

PV (N) -PP (N) -BETA*W (N) **2 .0/ (AV**2 . 0*RHO) 

PVTTL(N)-PV(N) * (1+ (GAMMA-1) /2*XM2) **GAMM1 

IF(XM2.GT.0.99)  THEN 
WRITE (6, 70)N,XM2 

70  FORMAT(3X, 'N-',I4,3X, 'XM2-',F6.3) 

ENDIF 


165 


30  TX»TVREAL(N-1) 


DO  80  NX»1,100 

T2=-5567.0/ALOG10 (PV(N) *TX**0 . 5/2 . 29E11 ) 

IF(ABS(T2-TX) .LT.1.0)  GO  TO  90 

TX=T2 

IF (NX. EQ. 100)  WRITE(6,*)  'NX-lOO* 

80  CONTINUE 

90  TVREAL(N)=T2 

DENV (N) =PV (N) / (R*TVREAL (N) ) 

C 

CC  WRITE (6, 20) Z/PL, P2/P0,T,W^,XM22 
CC  20  FORMAT (IX, 5 (3X, F15. 10) ) 

5  CONTINUE 
C 

DO  50  N=1,KBURN 
TVREAL (N) -TVREAL (KBURN+1) 
PV(N)=PV(KBURN+1) 

PVTTL (N) =PVTTL (KBURN+1 ) 

50  CONTINUE 
RETURN 
END 


C 

C 

C 

C*  SUBROUTINE  TES  * 


SUBROUTINE  TES 


C  CKL  -THERMAL  CONDUCTIVITY  OF  LIQUID  LIH  (W/M-K) 

C  CKSS  -THERMAL  CONDUCTIVITY  OF  SOLID  LIH  (W/M-K) 

C  CKP  -NODAL  THERMAL  CONDUCTIVITY  (W/M-K) 

C  CL  -HEAT  CAPACITY  OF  LIQUID  LIH  (J/K-KG) 

C  CS  -HEAT  CAPACITY  OF  SOLID  LIH  (J/K-’"') 

C  CP  -NODAL  HEAT  CAPACITY  (J/K-KG) 

C  DENL  -DENSITY  OF  LIQUID  LIH  (KG/M3) 

C  DENS  -DENSITY  OF  SOLID  LIH  (KG/M3) 

C  DENP  -NODAL  DENSITY  (KG/M3) 

C  DR  -SPACE  INCREMENT  IN  RADIAL  DIRECTION  IN  TES  (M) 

r  DT  -TIME  STEP  (SEC) 

C  DTHETA  -ANGULAR  INCREMENT  (RTJIIANS) 

C  DZ  -SPACE  INCREMENT  IN  AXIAL  DIRECTION  (M) 

C  HF  ’SE  -FUSION  HEAT  OF  LIH  (J/KG) 
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c 

RO 

-OUTSIDE  RADIUS  OF 

HEAT  PIPE  WALL(M) 

c 

c 

R 

TMELT 

-NODAL  RADIUS  (M) 

-MELT  TEMPERATURE  Ot  LIH  (K) 

c 

T(I,K) 

-NODAL  TEMPERATURE 

(K) 

c 

TV(K) 

-VAPOR  TEMPERATURE 

(K) 

C%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  ■!.%“„% 
REAL  *8  T (40, 100) , TT(40, 100) ,TOLD (40, 100) ,TVOLD (100) , TV (100) 

REAL  *8  TSTAR(40, 100) ,TP,TW,TE,TS,TN,DTV(100) 

REAL  *8  DH(40, 100) ,FH(40,100) ,FT(40, 100) ,CSTAR(40, 100) 

REAL  *8  P(100),Q(100),CKTES1(100),F,F1,T1,T2 

COMMON  /HPIPE/  TV,DQDATE (100) , 2L,DZ,NZ,NDT,QVTES (100) , 

$  NTESTU , TIME , KBURN, NWRI TE , XWRI TE , NGROV , RW , NADIA, P V ( 1 0  0 ) 

COMMON  /PCM/  TVOLD,T,TT,DTV,DT,NADI,NR,DTHETA,RO, 

S  TSTART, NTHETA, QPCM, QMELT, QPCMTTL 


F=  0.01 
Fl=2.0-F 

IF (NADI. EQ. 4)  GO  TO  400 
IF (NADI. EQ. 6)  GO  TO  600 


C  TES  GEOMETRY 
DR=RO/NR 


C  THERMAL  CONDUCTIVITY,  HEAT  CAPACITY  AND  DENSITY  FOR  TES 


CKL=2. 1 
CKSS=4.2 

CKM=2 .0* (CKSS*CKL) / (CKSS+CKL) 

CL=73'70.0 

CS=6280.0 

CM= (CL+CS) /2 . 0 

DENL-550 . 0 

DENS=DENL 

DENM= (DENL+DENS) /2 . 0 
CSTEEL-500 . 0 


C  MFLTING  TEMPERATURE  AND  HEAT  FUSION 

TMELT=956 . 0 
HFU5E=2580000 . 0 

QPCMTTL=3 . 1416*R0**2. 0*ZL*DENM*HFUSE*NTESTU 

Tl=951.0 
T2=961 . 0 
H1=CS*T1 

H2=CS*TMELT+CL* (T2-TMELT) +HFUSE 
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C  TES  STARTING  TEMPERATURE 

IF(NDT.EQ.l)  THEN 
DO  30  K=1,NZ 
DO  36  1=1, NR 
T (I, K) =TSTART 
TOLD (I, K) =TSTART 
36  CONTINUE 
30  CONTINUE 
END  IF 


DO  50  1=1, NR 
DO  55  K=1,NZ 

R=RO- (I-O . 5) *DR 


C  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I 
TP=TOLD (I, K) 

IF (I . EQ. 1)  TN=TP+2* (TVOLD (K) -TP) • ( 
IF(I.EQ.NR)  TS=TP 
IF(K.EQ.l)  TW=TP 
IF(K.EQ.NZ)  TE=TP 


IF(I.NE.l)  TN=TOLD(I-l,K) 
IF(I.NE.NR)  TS=TOLD(I+l,K) 
IF(K.NE.l)  TW=TOLD (I,K-1) 
IF(K.NE.NZ)  TE=TOLD (I,K+1) 


IF(TP.LE.TMELT)  THEN 
IF(TN.LE.TMELT)  CKN=CKSS 
IF<TS.LE.TMELT)  CKS=CKSS 
IF(TW.LE.TMELT)  CKW=CKSS 
IF(TE.LE.TMELT)  CKE=CKSS 


IF (TN.GT.TMELT) 
IF(TS.GT.TMELT) 
IF(TW.GT.TMELT) 
IF (TE.GT.TMELT) 
END  IF 


CKN=CIU^ 

CKS=CKM 

CKW=CKM 

CKE=CKM 


IF(TP.GT.TMELT)  THEN 
IF (TN. LE.TMELT)  CKN=CKM 
IF(TS.LE.TMELT)  CKS=CKM 
IF (TW. LE.TMELT)  CKW=CKM 
IF (TE. LE.TMELT)  CKE=CKM 


IF(TN.GT.TMELT)  CKN=CKL 


1-DR/ (4*RO) 
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IF(TS.GT.TMELT)  CKS=CKL 
IF (TW.GT.TMELT)  CKW=CKL 
IF(TE.GT.TMELT)  CKE-CKL 
ENDIF 

DH (I,K) = (DT/ (DENM*R*DR**2. 0) ) * (CRN* (R+DR/2) * (TN-TP) 

$  +CKS* (R-DR/2) * (TS-TP) ) 

S  + (DT/ (DENM*DZ**2 . 0) ) * (CKW* (TW-TP) +CKE* (TE-TP) ) 


C  MELTING 

IF(DH(I,K) .GE.0.0)  THEN 
IF(TP.LT.Tl)  FH(I,K)=CS*TP 

IF  (TP . GE .  T1 .  AND .  TP  .  LE .  T2)  FH  (I ,  K)  »  (H2-H1 )  *  (TP-Tl )  /  (T2-T1 )  +H1 
IF(TP.GT.T2)  FH(I,K)=(TP-T2) *CL+H2 


H»FH(I,K)+DH(I,K) 


IF(H.LT.Hl)  TSTAR(I,K)«'H/CS 

IF  (H .  GE .  HI .  AND  .  H .  LE  .  H2 )  TSTAR  (I ,  K)  •=  (H-Hl )  *  (T2-T1 )  /  (H2-H1 )  +T1 
IF(H.GT.H2)  TSTAR(I,K)=(H-H2) /CL+T2 


IF(FH(I,K) .LT.Hl.AND.H.LT.Hl)  CSTAR ( I , K) =CS 

IF(FH(I,K) .GE.H1.AND.FH(I,K) . LE . H2 . AND . H . GE . HI . AND . H . LE . H2 ) 

$  CSTAR(I,K)«(H2-H1)/(T2-T1) 

IF(FH(I,K) .GT.H2.AND.H.GT.H2)  CSTAR (I , K) -CL 

IF( (FH(I,K) .LT.Hl.AND.H.GT.Hl) .OR. (FH(I,K) . LT . H2 . AND . H . GT . H2 ) ) 
S  CSTAR(I,K)-DH (I,K) / (TSTAR (I, K) -TP) 

ENDIF 


C  SOLIDIFICATION 

IF(DH(I,K) .LT.0.0)  THEN 

IF(TP.LT.Tl)  FH(I,K)=CS’*TP 

IF(TP.GE.T1.AND.TP.LE.T2)  FH ( I , K) = (H2-H1 ) * (TP-T2) / (T2-T1)+H2 
IF(TP.GT.T2)  FH(I,K)-(TP-T2) *CL+H2 


H-FH(I,K)+DH(I,K) 


IF(H.LT.Hl)  TSTAR(I,K)-H/CS 

IF(H.GE.H1.AND.H.LE.H2)  TSTAR  ( I ,  K)  »  (H-H2 )  *  (T2-T1)  /  (H2-H1)  1-T2 
IF(H.GT.H2)  TSTAR(I,K)-(H-H2) /CL+T2 
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IF(FH(I,K) .LT.Hl .AND.H.LT.Hl)  CSTAR (I, K) =CS 

IF(FH(I,K) .GE.H1.AND.FH(I,K) . LE . H2 . AND . H . GE . HI . AND . H . LE . H2 ) 

$  CSTAR(I,K)=(H2-H1)/(T2-T1) 

IF{FH(I,K) .GT.H2.AND.H.GT.H2)  CSTAR ( I, K) =CL 

IF( (FH(I,K) .GT. HI. AND.H.LT.Hl) .OR. (FH(I,K) . GT . H2 . AND . H . LT . H2 ) ) 
$  CSTAR ( 1 , K) =DH ( I , K) / (TSTAR ( I , K) -TP ) 

ENDIF 


55  CONTINUE 
50  CONTINUE 


C& &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&£&&&&&& 
C  &&&&&&&  COMPUTATION  FOR  THE  FIRST  HALF  TIME  STEP  &&&&&&&&& 
QPCM-0 . 0 
QMELT=0 . 0 


DO  180  1=1, NR 
DO  160  K-1,NZ 

R=RO-(I-0.5)*DR 

C  t  I  I  1  I  I  I  i  I  I  I  I  I 

TP=TOLD(l,K) 

IF(I .EQ. 1)  TN-TP+2* (TVOLD(K) -TP) » (1-DR/ (4*RO) ) 
IF(I.EQ.NR)  TS-TP 
IF(K.EQ.l)  TW-TP 
IF(K.EQ.NZ)  TE=TP 

IFd.NE.l)  TN-TOLD(I-l,K) 

IF(I.NE.NR)  TS-TOLD (I+1,K) 

IF(K.NE.l)  TW-TOLD(I,K-l) 

IF(K.NE.NZ)  TE=TOLD(I,K+l) 

IF(TP.LE.TMELT)  THEN 
IF(TN.LE.TMELT)  CKN=CKSS 
IF(TS.LE.TMELT)  CKS=CKSS 
IF(TW.LE.TMELT)  CKW=CKSS 
IF(TE.LE.TMELT)  CKE=CKSS 

IF(TN.GT.TMELT)  CKN-CKM 
IF(TS.GT.TMELT)  CKS-CKM 
IF(TW.GT.TMELT)  CKW-CKM 
IF (TE.GT.TMELT)  CKE=CKM 
ENDIF 

IF(TP.GT.TMELT)  THEN 
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IF(TN.LE.TMELT)  CKN-CKM 
IF(TS.LE.TMELT)  CKS-CKM 
IF(TW.LE.TMELT)  CKW-CKM 
I F { TE . LE . TMELT )  CKE-CKM 


IF(TN.GT.TMELT) 
IF (TS.GT. TMELT) 
IF (TW.GT. TMELT) 
IF (TE.GT. TMELT) 
ENDIF 


CKN-CKL 

CKS-CKL 

CKW-CKL 

CKE=CKL 


DENP-DENM 

CP-CSTAR(I,K) 


IF(I.EQ.l)  CKTESl (K)-CKN 


VOLP- ( <R+DR/2) **2 . 0- (R-DR/2) **2.0) *3. 1416*DZ*NTESTU 

IF(T(I,K) .LT.Tl)  THEN 
QPCM»QPCM+DENS‘VOLP*CS* (T (I, K) -TSTART) 

QMELT-QMELT 

ENDIF 

IF(T(I,K) .GE.T1.AND.T{I,K) .LE.T2)  THEN 
QPCM-QPCM+DENS*VOLP*CS* (T1 -TSTART) +DENM*VOLP* (FH (I, K) -HI) 
QMELT«QMELT+DENM*VOLP* (FH(I,K) -HI) 

ENDIF 

IF(T{I,K) .GT.T2)  THEN 

QPCM-QPCM+  DENS*VOLP*CS* (TMELT-TSTART) +DENM*VOLP*HFUSE 
$  +DENL*VOLP*CL* (T(I,K) -TMELT) 

QMELT-QMELT+DENM*VOLP*HFUSE 
ENDIF 


Cl  I  I  I  I  I  I  I  I  I  I  I  I  I 
C— CKW*F1 

A-2 . 0*DENP*CP*DZ**2 . 0/DT+ (CKW+CKE) *F1 
B— CKE*F1 

D2- (DZ**2.0*CKN* (R+DR/2.0) /R/DR**2.0) *F 
D3- (DZ**2.0*CKS* (R-DR/2. 0) /R/DR**2. 0) *F 
Dl-2 . 0*DENP*CP*DZ**2 . 0/DT-D2-D3 

IF(I .EQ. 1)  D8- (2 . 0*DZ**2.0*CKN* (R+DR/4) /R/DR**2 . 0) *F 


IF(I.EQ.l.OR.I.EQ.NR)  THEN 

IFd.EQ.l  )  D-(D1+D2-D8)*T(I,K)+D3*T(I  +  1,K)+D8*TV0LD(K) 
IF(I.EQ.NR)  D-(D1+D3)*T(I,K)+D2*T(1-1,K) 

ELSE 

D-D1*T(I,K)+D2*T (I-1,K)+D3*T(I+1,K) 

ENDIF 


171 


IF(Dl.LT.O.O)  WRITE (6, 128) I, K 
128  FORMAT (IX, 'Dl-(-)  IN  ROOP3  I, K-' , 2 (2X, 13) ) 
IF(Dl.LT.O.O)  GO  TO  1000 

IF(K.EQ.l)  THEN 
C=0.0 

A-A-CKW*F1 

B-B 

D-D 

ENDIF 

IF(K.EQ.NZ)  THEN 
C-C 

A-A-CKE*F1 

B-0.0 

D-D 

ENDIF 

C 

C  TDMA  ALGORITHM 
C 

IF(K.EQ.l)  THEN 
P(K)— B/A 
Q(K)-D/A 
ELSE 

P(K)— B/  (A+C*P(K-1)  ) 

Q(K)-(D-C*Q(K-1)  )  /  (A+C*P{K-1) ) 

ENDIF 


IF(NDT.EQ.l)  THEN 

QVTES (K) -NTESTU*CKTES1 <K) *(2.0*3. 1416* (RO-DR/4 ) *DZ) 
$  *(T(I,K)-TVOLD(K))/(DR/2) 

DQDATE(K) -QVTES (K) / (NTESTU*2 . 0*3 . 1416*RO*DZ) 

ENDIF 


160  CONTINUE 


DO  170  K1-1,NZ 
K-NZ+l-Kl 

IF(K.EQ.NZ)  TT(I,K)=Q(K) 

TF(K.NE.NZ)  TT(I,K)=P (K) *TT ( I , K+1 ) +Q (K) 
170  CONTINUE 
180  CONTINUE 

DO  193  K-1,NZ 
DO  195  I-1,NR 
T(I,K)-TT(I,K) 

195  CONTINUE 
193  CONTINUE 
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IF (NADI. EQ. 123)  GO  TO  1000 
C 

££££££££  COMPUTATION  FOR  THE  SECOND  HALF  TIME  STEP  ££££££££ 

400  DO  380  K-1,NZ 
DO  360  I-1,NR 

R-RO-(I-0.5)*DR 

C  I  I  I  I  I  I  I  I  I  I  I  I  I 

TP-TOLD (I, K) 

IFd.EQ.  1)  TN-TP+2*  (TVOLD(K)  -TP)  *  (1-DR/  (4*RO)  ) 

IF(l.EQ.NR)  TS-TP 
IF(K.EQ.l)  TW-TP 
IF(K.EQ.NZ)  TE-TP 

IF(I.NE.l)  TN-TOLD(I-l,K) 

IF(I.NE.NR)  TS-TOLD(I+l,K) 

IF(K.NE.l)  TW-TOLD(I,K-l) 

IF(K,NE.NZ)  TE-TOLD(I,K+l) 

IF(TP.LE.TMELT)  THEN 
IF(TN.LE.TMELT)  CKN-CKSS 
IF(TS.LE.TMELT)  CKS-CKSS 
IF(TW.LE.TMELT)  CKW-CKSS 
IF(TE.LE.TMELT)  CKE-CKSS 

IF(TN.GT.TMELT)  CKN-CKM 
IF(TS.GT.TMELT)  CKS-CKM 
IF (TW.GT.TMELT)  CKW-CKM 
IF(TE.GT.TMELT)  CKE-CKM 
ENDIF 

IF(TP.GT.TMELT)  THEN 
IF(TN.LE.TMELT)  CKN-CKM 
IF(TS.LE.TMELT)  CKS-CKM 
IF(TW.LE.TMELT)  CKW-CKM 
IF(TE.LE.TMELT)  CKE-CKM 

IF(TN.GT.TMELT)  CKN-CKL 
IF(TS.GT.TMELT)  CKS-CKL 
IF (TW.GT.TMELT)  CKW-CKL 
IF (TE.GT.TMELT)  CKE-CKL 
ENDIF 


DENP-DENM 
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CP-CSTAR(I,K) 


Cl  I  I  I  I  I  I  I  I  I  I  I  I  I 

C=-CKN* (R+0 . 5*DR)  *F1 

A»2 . 0*DENP*CP*R*DR**2 . 0/DT+ (CKN* (R+0 . 5*DR) 
$+CKS* (R-0 . 5*DR) ) *F1 
B— CKS*  (R-0. 5*DR)  *F1 
D4- (R*DR**2 . 0*CKW/D2**2 . 0) *F 
D5- (R*DR**2 . 0*CKE/D2**2 . 0) *F 
IFCK.EQ.l)  D4-0.0 
IF(K.EQ.N2)  D5=0.0 
Dl=2 . 0*DENP*CP*R*DR**2. 0/DT-D4-D5 

IF(K.EQ.1.0R.K.EQ.N2)  THEN 
IF(K.EO.l)  D-D1*T(1,K)+D5*T(I,K+1) 
IF(K.EQ.N2)  D=D1*T(I,K)+D4*T(I,K-1) 

ELSE 

D=D1*T(I,K)+D4*T(I,K-1)+D5*T{I,K+1) 

ENDIF 

IF(Dl.LT.O.O)  WRITE <6, 328) I, K 
328  FORMAT (IX, 'Dl=(-)  IN  ROOPl  I, K-' , 2 (2X, 13) ) 
IF(Dl.LT.O.O)  GO  TO  1000 

IF(I.EQ.l)  THEN 
C-0.0 

A-A+CKN*R*F1 

B-B 

D=D+2 . 0*CKN* (R+DR/4) *TV (K) *F1 
C  D-D+2 . 0*CKN* (R+DR/4) * (TVOLD (K) +DTV (K) ) *F1 

ENDIF 

IF(I.EQ.NR)  THEN 
C-C 

A-A-CKS* (R-DR/2 . 0 ) *F1 
B-0.0 
D-D 
ENDIF 
C 

C  TDMA  ALGORITHM 
C 

IFd.EQ.l)  THEN 
P(I)— B/A 
Q(I)-D/A 
ELSE 

P(I)— B/  (A+C*P(I-1)  ) 

Q(I)»=(D-C*Q(I-1)  )  /  (A+C*P(I-1)  ) 

ENDIF 


360  CONTINUE 
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DO  370  11=1, NR 
I-NR+l-Il 

IF(I.EQ.NR)  TT(I,K)=Q(I) 

IF(I.NE.NR)  TT(I,K)-P (I) *TT (I+l , K) +Q (I) 
370  CONTINUE 
380  CONTINUE 


DO  310  K-1,NZ 

QVTES (K) -NTESTU»CKTES1 (K) * (2 . 0* 3 . 1416* {RO-DR/4) *DZ) * 
$  <TT(1,K)-TV{K))/<DR/2.0) 

310  CONTINUE 

IF(NADI.EQ.4)  GO  TO  1000 


600  DO  393  K-1,NZ 
DO  395  I-1,NR 
T(I,K)-TT(I,K) 
TOIiD(I,K)-TT(I,K) 
395  CONTINUE 
393  CONTINUE 


1000  RETURN 
END 


C 

Q  ********************** 

C  *  SUBROUTINE  HWANGBO  * 
^  ********************** 


SUBROUTINE  HWANGBO 


C  ONE  DIMENSIONAL  CAPILLARY  LIQUID  FLOW 


C 

C 

C 

c 

c 


QVL(K)  -  EVAPORATION  (-)  OR  CONDENSATION  (+)  RATE  PER  UNIT  HEAT 
PIPE  LENGTH  (W/M) 

VISLIQ  -  LIQUID  DYNAMIC  VISCOSITY  (NS/M2) 

TENLIQ  -  LIQUID  SURFACE  TENSION  (N/M) 

LAMBDA  -  LIQUID  LATENT  HEAT  (J/KG) 


REAL  *8  TV(IOO) 


DIMENSION  MDOT(IOO) 

DIMENSION  ALPHA(lOl) ,ALPHAD(101) ,OMEGA(101) , SAT(lOl) 
DIMENSION  WIDTH(lOl) , DEPTH (101) , AO (101) ,A1 (101) ,QTERM(101) 
DIMENSION  VOLFLW(lOl) ,DVDALP (101) ,W(5) ,V(5) ,QVL(100) 
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COMMON  /HPIPE/  TVjDQDATEEdOO)  ,  ZL,DZ,NZ,NDT,QVTES  (100)  , 
$NTESTU, TIME, KBURN, NWRITE, XWRITE, NGROV, RW, NADIA, PV (100 ) 

COMMON  /LIQ/  QV(IOO) ,PL(100) ,SUMQVI 
$  CL, NZPl , WIDBAR, DEPTHG,  DENL, VISLIQ, TENLIQ, LAMBDA 


Cl  KBURN- 0 

DO  3  K-1,N2 

Kl-NZ-K+1 

QVL(Kl)— 1.0*QV(K)  /NGROV/DZ 
3  CONTINUE 

C  ******  boundary  and  initial  CONDITIONS  ****** 

DO  5  K-1,NZP1 
DEPTH (K) -DEPTHG 
WIDTH (K) -WIDBAR 
5  CONTINUE 

C  GUESS  ALPHA (1) ,  OMEGA (1) 

ALPHA{1)=3. 1416/2.0 
OMEGA (1) -0.0 

DO  10  K=1,NZP1 
X-DEPTH (K) /WIDTH (K) /2 . 0 

AO (K)— 0.7058+0. 8685*X+0.1646*X**2. 0-0. 0145*X**3.0 
A1 (K) -0.1624+0. 3167*X-0. 1075*X**2. 0+0. 0073+X**3.0 
QTERM  <K) -VISLIQ/ (TENLIQ*DENL*WIDTH (K) **3.0) 

VOLFLW (K) — (AO (K) +A1 (K) *ALPHA (K) ) 

DVDALP(K)— A1  (K) 


IF(ALPHA(K) .LT.  (3.1416/180.0) )  THEN 
KBURN-NZPl-K 
GO  TO  30 
ENDIF 


IF(K.EQ.NZPl)  GO  TO  30 

VOLF-VOLFLW(K) 

DVDA-DVDALP (K) 

W(l)-0.0 

V(l)-0.0 

DO  20  N-1,4 

ALP-ALPHA (K) +0 . 5*DZ*W (N) 

OME-OMEGA (K) +0 . 5*DZ*V (N) 

W (N+1 ) -OMEGA (K) +0 . 5*DZ*V (N) 

V(N+1)  —  (DVDA/VOLF+1 .0/TAN(ALP)  )  *  (ABS  (OME)  )  **2. 0 
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$  +QTERM{K) * (QVL(K) /LAMBDA) / (VOLF*SIN (ALP ) ) 

20  CONTINUE 

ALPHA(K+l)-ALPHA(K)  +  (DZ/6.0)*  (W(2)+2.0’*W(3)+2.0*W(4)+W(5)  ) 
OMEGA (K+1) -OMEGA (K)+(DZ/6.0) * (V (2) +2 . 0*V (3) +2 . 0*V (4 ) +V (5) ) 


10  CONTINUE 


30  DO  50  K=1,NZP1-KBURN 

MDOT  (K)  -VOLFLW  (K)  *SIN  (ALPHA  (K)  >  *OMEGA  (K)  /QTERM  (K) 
50  CONTINUE 


IF( (NDT/NWRITE) .EQ. (NDT/XWRITE) )  THEN 
CC  WRITEE(6,*)'  • 

CC  WRITEE{6,*)' KBURN- ' , KBURN 


CC 

WRITES (6, *)  • 

t 

CC 

WRITEE(6, *)  ' 

K 

ALPHAD 

OMEGA 

MDOT 

DO  60  K-1,NZP1 -KBURN 
ALPHAD{K)- (180. 0/3. 1416) ‘ALPHA (K) 

CC  WRITEE(6, 65) K, ALPHAD (K) , OMEGA (K) , MDOT (K) 

CC65  FORMAT(1X,I4,2(2X,F10.5),2X,F11.6) 

60  CONTINUE 
ENDIF 
RETURN 
END 
C 


Q  A  ****  tk  A  A  ************  * 

C  *  SUBROUTINE  LIQUID  * 

Q  ********************* 

c 

c 

SUBROUTINE  LIQUID 

C  ONE  DIMENSIONAL  CAPILLARY  LIQUID  FLOW 

C  VISLIQ  -  LIQUID  DYNAMIC  VISCOSITY  (NS/M2) 

C  TENLIQ  -  LIQUID  SURFACE  TENSION  (N/M) 

C  LAMBDA  -  LIQUID  LATENT  HEAT  (J/KG) 

REAL  *8  TV(IOO) ,KPERM 

DIMENSION  QL(IOO) ,DELPL(100) ,DPVPL(100) 
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COMMON  /HPIPE/  TV,DQDATE(100) , ZL,DZ,NZ,NDT,QVTES (100) , 
$NTESTU, TIME, KBURN, NWRITE, XWRITE, NGROV, RW, NADIA, PV (100) 

COMMON  /LIQ/  QV(IOO)  ,  PLdOO)  ,  SUMQVI,KPVMIN,PCMAX,DPVPLMAX, 

$  CL,  NZPl , WIDBAR, DEPTHG, DENL, VISLIQ, TENLIQ, LAMBDA 


FLREL=14.3 

RM-RW+DEPTHG/2 

AW=2 . 0*3 . 1416*RM*DEPTHG 

EPS=NGROV*WIDBAR/ (2 . 0* 3 . 1416*RM) 

RHL=2 . 0»W1DBAR*DEPTHG/ (2 . 0*WIDBAR+2 . 0*DEPTHG) 
KPERM=2 . 0*EPS*RHL**2 . 0/FLREL 

FL=VISLIQ/ (KPERM*AW*DENL*LAMBDA) 


QL (KBURN+1 ) =QV (KBURN+1 ) 

DO  10  K=KBURN+2, NADIA 
QL (K) »QL (K-1 ) +QV (K) 

10  CONTINUE 


DO  20  K=  NZ,  KBURN+1,  -1 

IF  (K.GT. NADIA)  PL(K)-=PV(K) 
IF (K.LE. NADIA)  THEN 
DELPL (K) =-FL*QL (K) *DZ 
PL(K)=PL(K+1)+DELPL(K) 
ENDIF 

20  CONTINUE 


DO  30  K-KBURN+1, NADIA 
DPVPL (K) -PV (K) -PL (K) 

30  CONTINUE 

DPVPLMAX-DPVPL (KBURN+1 ) 


C  IF (DPVPL (KBURN+1 ) .GE.PCMAX)  THEN 

C  DO  60  K-KBURN+1, NADIA 

C  IF (DPVPL (K) .GT.PCMAX)  KBURNNEW-K-1 

C  60  CONTINUE 
C  ENDIF 

C  IF ( DPVP  L ( KBURN+ 1 )  . LT . PCMAX . AND . KBURN . EQ . 0 )  KBURNNEW-KBURN 

C  IF (DPVPL (KBURN+1 ) .LT. PCMAX. AND. KBURN. NE. 0)  KBURNNEW-KBURN- 1 

C 

C  KBURN-KBUBNNEW 
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KBURN=0 


RETURN 

END 
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